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ABSTRACT 

We investigate whether coagulation models of planet formation can explain 
the observed size distributions of transneptunian objects (TNOs). Analyzing 
published and new calculations, we demonstrate robust relations between the 
size of the largest object and the slope of the size distribution for sizes 0.1 km 
and larger. These relations yield clear, testable predictions for TNOs and other 
icy objects throughout the solar system. Applying our results to existing obser- 
vations, we show that a broad range of initial disk masses, planetesimal sizes, and 
fragmentation parameters can explain the data. Adding dynamical constraints 
on the initial semimajor axis of 'hot' KBOs along with probable TNO formation 
times of 10-700 Myr restricts the viable models to those with a massive disk 
composed of relatively small (1-10 km) planetesimals. 

Subject headings: Planetary systems - Planets and satellites: formation - Planets 
and satellites: physical evolution - Kuiper belt: general 



1. INTRODUCTION 

The transneptunian objects (TNOs, small icy planets orbiting the Sun beyond Neptune) 
provide crucial tests of planet formation theories. With its rich dynamical structure, the cur- 
rent orbital architecture of TNOs constrains models for the dynamical interactions among 
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newly-formed gas giants and for the origin of the Oort cloud (e.g., iMorbidelli et al.ll2008l ). 
The diverse populations o f binary TNOs provide additional cons t raints on the dynamical his- 
tory of the solar system ( iNoll et al.ll2008l ; iNesvorny et al.ll201ll ; iMurray-Clay fc Schlichting 
201ll ). Finally, the sizes of the largest TNOs and the size distributions within the dynam- 
ical families of TNOs constrain the formation histories of the gas giants, the TNOs, and 



perhaps the terrestria l planets (e.g., iKenyon fc Luulll998l ; lKenyonll2002l ; iGomes et al.l 12005 
Morbidelli et~alll2008h . 



There are two broad concepts for the formation and evolution of TNOs. Both begin 
with small dust grains within a massive gaseous disk surrounding a young star. As the disk 



evolves, dust grains collide, merge, and grow into mm- to cm-sized ob 
from the gas and fall into the midplane of the disk (e.g.. iBirnstiel et al- 



erts which decouple 



2010. and references 



therein). Once a large fraction of the solid material is in the midplane, the two paths to 
TNOs diverge. In one path, icy objects continue to collide, merge, and grow into roughly 
km-sized planetesimals and then into much larger TNOs (e.g., IWeidenschillingi ll 99 71 ) . In the 
other path, instabilities or turbulent eddies in the dis k concentrate icy objects into massive 



clumps which collapse directly into large TNOs (e.g., IJohansen et al.l 120071 ; IChambera 12010 



Cuzzi et al.l |2010| ; iPan et al.l 120111 ). Once TNOs form, the evolutionary paths converge: 
dynamical interactions with the gas giants p roduce the diverse architecture of the current 



TNO population (e.g.. Morbidelli et al.ll2008f ). 



Testing the two formation paths requires robust theoretical models which make clear 
predictions for the maximum radius and the size distribution of TNOs. At present, analytic 
theory and numerical s i mulations cann ot predict the outcome of instability mechanisms 
( iChiang fc Youdinll2010l : IPan et al.ll201ll ). However, coagulation models yield robust predic- 
tions for the time evo l ution of the size distribution and the sizes of the largest obje cts (e.g., 
Kenyon fc Luulll999al ; lKenyonll2002l ; IKenyon fc Bromley! 12004a ; IKenyon et al.ll2008l and ref- 
erences therein). Thus, it is possibl e to compare the resul t s of coagulation calculations with 
the o bservable properties of TNOs (IKenyon fc Luulll999bl ; IKenyon fc Bromleyll2004d ; iFraser 
2009h . 



To test the coagulation model in detail, we analyze a set of published ( IKenyon fc Bromley 



20081 . |2010| ) and new coagulation calculations. Our numerical models follow the collisional 
and dynamical evolution of small planetesimals into TNOs at orbital semimajor axes a 
= 15-150 AU for evolution times 1-10 Gyr. This r ange of semimajor axes lies outside 
the most likely region of giant planet formation (e .g. , iRafikovl l2004t iGoldreich et al.l 12004 ; 
Kennedy fc Kenyonll2008l ; IKenyon fc Bromley! 120091 ). By considering a broad range of initial 
planetesimal sizes (and size distributions) and a range of fragmentation parameters, we de- 
rive the time evolution of the radius of the largest object and the slope of the size distribution 
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as a function of initial conditions and distance from the Sun. Comparisons of our results 
with the observations allows an assessment of the match between the observed properties of 
TNOs and the models. 

Coagulation models yield robust predictions for the relation between the size of the 
largest object and the slope of the size distribution. Although this relation is insensitive to 
distance from the Sun and fairly insensitive to the fragmentation parameters, the evolution 
depends on the initial disk mass, the initial size of the largest planetesimals, and the initial 
size distribution of planetesimals. 

Using recent observations from large surveys of TNOs, we demonstrate that collisional 
evolution often yields size distributions similar to those observed. On their own, the observa- 
tions favor models with 1-10 km planetesimals relative to those with 100 km planetesimals. 
Adding dynamical constraints on the origin of the hot and cold populations of KBOs to our 
analysis, we show that models in (i) massive disks with 1 km planetesimals and (ii) low mass 
disks with 10 km planetesimals provide equally satisfactory matches to the data. Requiring 
that TNOs form prior to the Late Heavy Bombardment eliminates models with low mass 
disks. Thus, all of the available data taken together imply that TNOs formed in a massive 
disk composed of 1 km planetesimals. 

We begin our discussion with an overview of the numerical code in £j2j After describing 
basic results of the calculations in §21 we compare the results with observations in §U We 
conclude with a brief summary in §5J 



2. PLANET FORMATION CALCULATIONS 

To calculate the formation and evolution of TNOs, we use a hybrid multiannulus 
coagulation- A-body code. We compute the collisional evolution of an ensemble of planetesi- 
mals in a circumstellar disk orbiting a central star of mass M*. The code uses statistical algo- 
rithms to evolve the mass and velocity distributions of low mass objects with time and an N- 



body algorithm to follow the individual trajectorie s of massive objects. iKenvon fc Bromley 



(120081 ) describe the statistical (coagulation) code; iBromley fc Kenyonl (120061 ) describe the 



A-body code. Here, we briefly summarize the basic aspects of our approach. 

We perform calculations on a cylindrical grid with inner radius a in and outer radius a out 
surrounding a star with M* = 1 M . The model grid contains A concentric annuli with 
widths 5ai = 0.025aj centered at semimajor axes dj. Calculations begin with a cumulative 
mass distribution n c (mik) oc m ifc 9in4t of planetesimals with mass density p p = 1.5 g cm -3 
and maximum initial mass m . Although our code explicitly tracks the mass distribution 
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of solid objects, in the rest of the text we describe initial conditions and results in terms of 
cumulative size distributions, n c (r). For a power law cumulative size distribution, n c oc r~ q , 
the power-law exponent is a factor of 3 larger than the exponent for the mass distribution. 
Thus, q init = 3q' init . For this paper, we adopt q init = 0.5 (most of the mass in the largest 
planetesimals) and 3.0 (mass distributed equally per logarithmic interval in mass/radius). 

Planetesimals have horizontal and vertical velocities h ik (t) and v^t) relative to a 
circular orbit. The horizontal velocity is related to the orbital eccentricity, e 2 k (t) = 1.6 
{hik{t)/VK,i) 2 , where Vk,% is the circular orbital velocity in annulus i. The orbital inclination 
depends on the vertical velocity, i 2 k (t) = sm~ 2 (2(vik(t)/VK,i))- 

The mass and velocity distributions of the planetesimals evolve in t ime due to inelastic 
collisions, dra g forces, and gravitational encounters. As summarized in iKenyon &: Bromley 
(j2004al . 120081 ). we solve a coupled set of coagulation equations which treats the outcomes of 
mutual collisions between all particles with mass rrij in annuli a^. We adopt the particle- 
in-a-box algorithm, where the physical collision rate is navf g , n is the number density of 



objects, a is the geometric cross-section, v is 


the relative ve 


ocity, and f n is the gravitational 


focusing factor ( 


Safronov 


1963; 


Lissauei 


1987 


; Spaute et al. 


1991: 


Wetherill & Stewart 


1993; 


Weidenschilline: et al. 


19971: 


Kenvon & Luu ] 


L998|; 


Kenvon & Bromlev 


2008|). For a specific 



mass bin, our solutions include terms for (i) loss of mass from gas drag and mergers with 
other objects and (ii) gain of mass from collisional debris and mergers of smaller objects. 

Aside from deriving the cross-section and relative collision velocity, the main ingredient 
in this calculation is the gravitational focusing factor for a collision between two planetesimals 
with masses (radii) and rrij (r$ and rj). In the high velocity, dispersion-dominated regime, 
we adopt a variant of the piecewise analytic approximation of ISpaute et al.l (119911 . see also 
Kenyon & Luu 1998): 



fg,a 



1 + P(v e /v) 2 
42.4042 (ve/v) 1 - 2 
10651.453 



v > 0.32v e , 

0.01 < v/v e < 0.32, 

v/v e < 0.01, 



where is an order unity con stant which accounts for the distribution of impact velocities 
(IGreenzweig fc Lissauerll 19921 ). v is the velocity of a planetesimal, and v e is the escape velocity 
of the merged pair of planetesimals. 

For the low velocity, shear-dominated regime, we define the angular velocity Q = 
(GM./d 3 ) 1 / 2 , the Hill radius r H = a(m/3M ) 1 / 3 , the Hill velocity v H = Qr H , and an angular 
scale factor ip = (pe/pp) 1 ^ 3 R©/o, where p Q is the mean density of the Sun and R is the 
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solar radius. In the iGreenberg et al.l ( 119911 ) approach, 

f (l+v 2 Jv 2 T )(v T /v)(r H /r T ) 
l9 ' s \ 0.5(1 + v 2 /v 2 T y/ 2 (v T /v)(r H /r T )(H/( ri + r,-)) 



v<v H , v < ip 1/2 v H 



(2) 



where ry = a((m, + rrij) /m*) 2 ^ is the Tisserand radius, uy = l.lfWay i s the Tisserand 



veloci ty, and if is the vertical scale height of the planetesimals. In the iGoldreich et al. 
( 120041 ) approach, 



if> 1 (vh/v) 
a' 3 / 2 



< 



<<A 1/2 



(3) 



Tests using both approaches yield similar results for the collision rate. The simpler IGoldreich et al 
( 120041 ) relations are computationally cheaper. For either approach, we employ a simple algo- 



rithm to affect a smooth transition in f g between the dispersion and shear rates ( iKenyon fc Luu 



19981 ). Because ip decreases with a, icy planet formation at 15-150 AU is less likely to enter 
the shear-dominated regime than rocky planet formation at 1 AU. 

Collision outcomes depend on the ratio Q c /Q*d, where Q* D is the collision energy needed 
to eject half the mass of a pair of colliding planetesimals to infinity and Q c is the center of 
mass collision energy. For two colliding planetesimals with masses mi and m-i, the mass of 
the merged planetesimal is 

m = mi + m2 - m e j , (4) 
where the mass of debris ejected in a collision is 



0.5 (m 1 + m 2 ) [ 



9/8 



(5) 



To place the debris in our grid of mass bins, we set the mass of the largest c ollision frag- 
ment as mi = 0.2m esc and adopt a cumulative size distribution n c oc r~ qd (IDavis et al. 
19851 ) . Here, we adopt = 2.5, the value expected for a system in collisional equilib 



20031 ; iKobayashi fc Tanaka 



201 



rium ( Dohnanvi 1969; William sfc Wetherilllll994l ; lTanaka fc IdalliggGuO'Brien fc Greenberg 



This approach allows us to derive ejected masse s for catastrophic co 



Q* D and for cratering collisions with Q r <C Q* D (see alsolWetherill &: Stewart 



1994 : 

200a 

e.g., 



Tanaka Ida 



1996 



Kobavashi &: Tanaka 



Stern fc Colwelllll997al ; lKenyon fc Luu 



1999a 



isions with Q r 



1993 



Williams fc Wetherill 



O'Brien fc Greenberg 



201 0j). Consistent with N-body simulations of collis ion outcomes 



Benz fc Asphaugj|l999l ; iLeinhardt et al.ll2008l ; iLeinhardt Stewartl 120091 ) , we set 



Q* D = Q b r A + Q gP y 



(6) 



where Qbr l3b is the bulk component of the binding energy, Q g p g T^ 9 is the gravity component 
of the binding energy, and r is the radius of a planetesimal. 
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To explore the sensitivity of our r esults to the fragmentation algorithm, we consider three 
sets of parameters /j (Fig. [Q. As in iKenyon fc Bromley! (120081 . 2010), strong planetesimals 



have L = (Qh = 10 1 , 10 3 , or 10 5 erg g" 1 , (3 b = 0, Q fi 



= 2.25 erg g 
M = 2 x 10 5 



2 cm l-75 



A, 



1.25; 



erg g 



- 1 cm - 4 , P b 



Benz fc Asphauglll999l ); weaker planetes imals have f„ 
—0.4, Q g = 0.33 erg g -2 cm 1 - 7 , f3 g = 1.3; iLeinhardt fc Steward 120091 ) . Here, we also consider 
a set of very strong planetesimals w ith f vs = (Qp = 1 1 , 10 3 , or 10 5 erg g -1 , /3& = 0, Q g 



10 4 erg g 2 cm, 



f3 g = 2.0; e.g., iDavis et al.l Il985h . At small sizes, these parameters 



are broadly consiste nt with results from laboratory experiments of impacts with icy targets 
and p rojectiles (e.g. 
20051 ). 



Rvan et al.lll999l : lArakawa et al.ll2002l ; iGiblin et al.ll2004J : iBurchell et al. 



For our calculations, these parameters have several consequences. When fragmentation 
begins, most of the mass is in planetesimals with r p=j 1-100 km. With Q* D (r = 1-100 km) > 
Q* D (r = 1-10 5 cm), the gravity regime of the fragmentation law sets the collision velocity for 
the onset of fragmentation. To achieve catastrophic fragmentation with Q c /Q*d = 1> weak 
planetesimals require the smallest collision velocity; very strong planetesimals require the 
largest. In all cases, once large objects begin to fragment, smaller objects also fragment. 
Thus, the onset of fragmentation for 1-100 km planetesimals initiates a collisional cascade 
where all small objects are ground to dust. 

To compute the evolution of the velocity distribution, we include collisional damping 
from inelastic collisions, gas drag, and gravitational interactions. For inelastic and elas- 



tic collisions, we fol low the statistical, Fokker-Planck approaches of lOhtsukil (119921 ) and 
Ohtsuki et al.l (120021 ). which treat pairwise interaction s (e.g., dynamical f r iction and vis- 
cous stirring) between all objects in all annuli. As in IKenyon fc Bromley! (120011 ). we add 
terms to tr eat the probability that ob jects in annulus i\ interact with objects in annulus 
12 (see also IKenyon fc Bromlev! l2004bl. KB08). We also compute long-range stirring from 
distant oligarchs (jWeidenschillingl Il989[ ) . For gas drag, we derive velocit y damping and the 



inward drift of material in the quadra tic, Stokes, and Epstein limits (e.g.. lAdachi et al.lll976 



Weidenschillindll977d : lRafikovil2004f ) 



To evolve the gas in time, we consider a si mple nebular model for the gas density. We 



adopt a scale height H, 



gas\ 



H, 



gasfl 



[a/a 



,1.125 



(IKenyon fc Hartmannlll987l ) and assume the 



gas surface density declines exponentially with time 



->gas 



[a,t) 



-'gas 



-£m O 



t /tgas 



(7) 



where S sas>0 and x m are scaling factors and t gas = 10 Myr is the gas depletion time. Although 
longer than the 2-5 Myr timescale estimated from ob s ervations of the lifetimes o f accretion 
disks in pre-main sequence s tars (jCurrie et al.l |2009| ; Kennedy &: Kenyonl |2009| ; iMamajek 
20091 ; IWilliams fc Ciezall201ll ). shorter gas depletion times have little impact on our results. 
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In the iV-body code, we directly integrate the orbits of objects with masses larger than 
a pre-set 'promotion mass' m pro . The calculations allow for mergers among the iV-bodies. 
Additional algorithms treat mass accretion from the coagulation grid and mutual gravita- 
tional stirring of n-bodies and mass batches in the coagulation grid. For the applications in 
this paper, the few large objects capable of promotion into the iV-body code never contain a 
significant fraction of the mass in an annulus and never contribute significantly to the local 
stirring. To treat situations where a few large objects might impact the evolution, we set 
10 26 g. However, our calculations never produced more than a few iV-bodies which 



m 



pro 



remain on nearly circular orbits throughout their evolution. 

The initial c onditions for these calcu lations are appropriate for a disk with an age of 
< 1-2 Myr (e.g., IWilliams &: Ciezal l201ll . and references therein). We consider systems of 



N = 64 annuli in disks where the initial surface density of solid material follows a power law 
in semimajor axis, 



where ai is the central radius of the annulus in AU, n = 1 or 3/2, and oc yyi, is a scaling 
factor. For a standard gas to dust ratio of 100:1, E 9aSj0 = 100 E^^M*) . To explore a range 
of disk masses similar to the observed range among the youngest stars, we consider E^o 
= 30 g cm -2 and x m = 0.01-3. Disks with x m ~ 0.1 have masses similar to the median 



disk masses observed ar ound young stars in nearby dark clouds ([Andrews fc Williamsl 12005 
Williams fc Ciezal 120111 ) . Somewhat larger scale factors, x w . ~ 1 T corr espond to mode ls of 



the minimum mass solar nebula models of IWeidenschillingl (jl977bl ) and lHayashil (119811 ). 



Table [T] summarizes the model grids. For each set of x m , we choose the extent of the 
disk, a = 15-75 AU or 30-150 AU; initial radius of the largest planetesimal, tq = 1, 10, 
100 km; the initial eccentricity Cq = 10 -4 and inclination i = 5 x 10 -5 ; the power law 
exponent of the initial size distribution, q init = 3 (equal mass per log interval in mass) or 
Qinit — 0.5 (most of the mass in the largest objects); the fragmentation parameters, / s , f vs , 
or f w ; and the evolution time t max . To understand the possible ra nge of outcomes, we r epeat 



calculations 5-15 times for each combination of initial conditions. iKenyon fc Bromley! (12008 



2010) discuss some aspects of these calculations in the context of observations of debris disks. 
Our calculations follow the time evolution of the mass and velocity distributions of 



objects with a range of radii, rj 



r min to Tij 



The upper limit r max is always larger 



than the largest object in each annulus. To save computer time, we do not consider small 
objects which do not significantly affect the dynamics and growth of larger objects, r min = 
100 cm. Erosive collisions produce objects with r^- < r min which are 'lost' to the model grid. 
Lost objects are more likely to be ground down into smaller objects than to collide with 
larger objects in the grid (see lKenyon fc Bromley! 120021 . 2002b, 2004a). 
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EVOLUTION OF THE SIZE DISTRIBUTION 



At the start of each calculation, dynamical processes drive the evolution. Gas drag 
damps the velocities of the smallest objects. Dynamical friction damps the orbital velocities 
of the la rgest objects; dynamic a l friction and viscous stirr ing excite velocities of the smallest 
objects (IGoldreich et al.l 12004] ; iKenyon fc Bromley! 12008k 2010). With stirring timescales 
much faster than growth timescales, it takes 10-50 orbits to achieve a relaxed distribution of 
velocities among all planetesimal sizes. For calculations with r$ = 1-10 km, this evolution 
occurs in the dispersion-dominated regime; dynamics and growth are slow and orderly. When 
r = 100 km, dynamical evolution rapidly drives the system from the shear regime to the 
dispersion regime. Thus, initial growth rates for large planetesimals are also slow and orderly. 

As the evolution proceeds, dynamical processes produce substantial gravitational focus- 
ing factors. Runaway growth begins. The largest objects then grow much more rapidly than, 
and 'run away' from, somewhat smaller protoplanets. The initial stages of this process are in 
the dispersion regime; accretion in the shear regime dominates the later stages. Throughout 
this evolution, the largest objects - with radii of 300-1000 km - continue to stir the smallest 
objects. Stirring reduces gravitational focusing factors, ending shear-dominated accretion 
and runaway growth. Oligarchic growth in the dispersion regime begins. 

Throughout oligarchic growth, protoplanets evolve in two distinct ways. The largest 
protoplanets - known as oligarchs - accrete material primarily from the reservoir of small 
planetesimals. Because gravitational focusing factors are small, growth is orderly. Small 
oligar chs grow faster than large oligarchs. Both grow faster than smaller planetesimals 
(e.g., IGoldreich et al.l 120041 ) . As all of the oligarchs grow, they stir the smallest planetes- 
imals to larger and larger orbital velocities. Eventually, collisions among small planetes- 
imals become destructive, producing smaller planetesimals and some debris. Collisions 
among smaller planetesimals are even more destru ctive, leading to a collisional cascade 
wher e small planetesimals are g r ound to dust (e.g., iDohnanvil Il969t IWilliams fc Wetherill 
19941 ; iKobayashi fc Tanakal |2010| ; iBelyaev fc Rafikovl 1201 ll . and references therein) . By rob- 
bing the oligarchs of small planetesimals, the collisional cascade limits the growth of the 



largest planets to ~ 3 - 10 x 10 3 km fe.g-. llnaba et al.ll2003l : Kenvon fc Bromlevll2008l . l2010l ; 



Kobavashi fc Tanakal l201Q[ ). 



On 1-10 Gyr tim e scales, the collisional cascade r emoves most of the i nitial mass in 
solids at a < 100 AU fcenvon & Bromlevl I2008I . boioh . iKenvon k Bromlevl fepQsh express 
the fraction of solid material remaining in the disk as f s « f Q Xm^(a/30 AU) 5 / 4 . Collisional 
growth proceeds more rapidly at smaller a in more massive disks. For a fixed evolution 
time, the cascade removes more mass at smaller a and less mass in less massive disks. 
For calculations with r = 1 km and strong (weak) planetesimals, /o,« ~ 0.12 (0.08) at 
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1 Gyr; / 0>s « 0.10 (0.07) at 10 Gyr flKenvon fc Bromley! 12008 . l2010h . For the very strong 
planetesimals considered here, fo )VS ~ 0.26 (0.16) at 1 Gyr (10 Gyr). Thus, the cascade 
removes 74% to 93% of the initial mass at 30 AU on 1-10 Gyr timescales. 

The amount of solid material remaining in the disk depends on the initial maximum size 
of the planetesimals. In calculations with ro = 10 km, the cascade removes roughly a factor 
of two less mass over 1-10 Gyr. When r = 100 km, the cascade typically removes a facto r 
of 7-8 less mass. Effective cascades require small planetesimals (IKenyon fc Bromley! 120101 ) . 
Ineffective cascades leave more mass in the disk and allow more massive planets to form. 
Thus, th e most massive planets grow in ensembles of large and very strong planetesimals 
(see also iKenvon fc Bromlev! 120081 boioh . 



Fig. [5] illustrates the outcomes of growth for an ensemble of 1 km planetesimals in eight 
annuli at 27-32 AU. Our starting conditions place most of the mass in 1 km planetesimals; 
the cumulative size distribution of planetesimals is nearly flat. During orderly growth, ~ 10 4 
planetesimals reach radii of nearly 10 km in 3 Myr. Once runaway growth begins, a handful 
of planetesimals grow as large as 100 km. In both cases, the size dist ribution follows a 



shallow power law at small radii and an exponential at large radii (e.g., iDavis et al.lll999 



Kenyon fc Luul Il999bl ). As the evolution makes the transition from runaway to oligarchic 
growth, the largest objects grow from ~ 100 km (10 Myr) to ~ 500 km (30 Myr) to ~ 2000 km 
(300 Myr). Because the oligarchs grow faster than smaller planetesimals, the size distribution 
becomes a fairly smooth power law from ~ 10 km to 
2004cl ; lFraseril2009l : ISchlichting fc Sarilboilh . 



2000 km (e.g., IKenyon &: Bromley 



At small siz es (r 5, 10 km), the collisional cascade produces inflection points in the 
size distribution (jKenvon fc Bromley! l2004q ; iPan fc Sarill2005l ; iBenavidez fc Campo Bagatin 



20091 ; iFraserl 120091 ) . Initially, the cascade 'erodes' planetesimals with r < 1 km. Smaller 



planetesimals erode more than large planetesimals. These planetesimals grow smaller with 
time. Because the cascade removes smaller planetesimals from an annulus more rapidly 
than larger planetesimals, the size distribution becomes shallower, creating an inflection 
point where erosion dominates growth. Eventually, larger collision velocities cause small 
planetesimals to shatter completely on impact. This process creates another inflection point 
where shattering dominates erosion. At the smallest sizes, debris from collisions of larger 
planetesimals dominates shattering, resulting in a sharp rise i n the size distribution (e.g., 
Kenvon fc Bromlev! lioolcl I2OO5I : IPan fc Sarilbood : lFraseril2009h . 



To analyze the time-evolution of the size distribution for a variety of initial conditions, 
we perform a series of least-squares fits to the results of our calculations (Fig. [3]). For 
each output size distribution at time U, we construct cumulative size distributions n c {r) for 
sets of eight annuli centered at distances a,- from the central star. With N = 64 annuli in 
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our calculations, we have eight cumulative size distributions at each timestep. We adopt a 
cumulative size distribution in the form of a power law n c oc r~ q and derive the best-fitting 
exponents for planetesimals with radii r = 0.1-1 km (qi), r = 1-10 km (g 2 ), T = 10-100 km 
(qs), and r > 100 km (g 4 ). To minimize fluctuations due to Poisson statistics, we restrict 
these fits to bins with n c (r) > 10. 

Although we could fi t the size distribution for r > 1 km to the product o f a power law 
and an exponential (e.g., IWetherilllll990l ; iTanaka fc Nakazawalll994l ; iLeei 120001 ). several fac- 
tors influence our choice of piece-wise fits to the overall size distribution. In our calculations, 
we set the initial maximum sizes of planetesismals at r$ = 1, 10, and 100 km. Measur- 
ing slopes between these sizes allows us to understand how our choices impact our results. 
Mergers and d estructive collisions often produce mult iple inflection points in the size dis- 
tribution (e.g.. iKenyon fc Bromley! 12004a ; lFraserll2009l ). Fitting several distinct power laws 



allows us to identify the inflection points easily and to compare power laws at sizes smaller 
and larger than the inflection points. Observation s of TNOs are often wel l -fit by double 



powe r laws with a 'break' at r ~ 10-100 km (e.g.. iFuentes &: Holmanl 120081 ; iFuentes et al. 



2009; Fraser fc Kavelaars 2009; Fraseret al. 2010 



Sheppard &: Tr u j illol 1 2 1 Oh . Our q 3 and g 4 



exponents allow us to compare calculations directly to observations in the outer solar system. 

In the next sections, we consider the time-evolution of the four power law exponents and 
the size of the largest object in an annulus r max . In every annulus, the growth of planetesimals 
follows a similar path from orderly growth to runaway growth to oligarchic growth and the 
collisional cascade. Thus, r max serves as a proxy for the evolution tim e t and allows us to 
compare re sults in different annuli at similar evolutionary states (e.g., IKenyon fc Bromley 
20081 . 120101 ) . We begin with a discussion of the evolution of the size distribution at large 
sizes (§3.1) and then describe the evolution at smaller sizes in subsequent subsections. 



3.1. The Size Distribution of Large Objects 

Fig. H] illustrates the time evolution of r max for the calculation described in Fig. [2J In 
each annulus, the largest objects first grow slowly; r max grows roughly linearly with time. 
Once runaway growth begins, large objects grow exponentially; r max reaches ~ 1000 km 
in 3 Myr at 15-18 AU, 20 Myr at 22-27 AU, 50 Myr at 33-40 AU, and 200 Myr at 48- 
58 AU. After planets achieve these radii, the collisional cascade removes small planetesimals 
more rapidly than planets accrete them. Thus, planets grow much more slowly. Despite the 
different timescales, planets in all annuli reach the same maximum size of roughly 3000 km. 

Fig. shows how g 4 evolves with time. When the first objects with r > 100 km form, 
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there are many 100 km objects and only a few with larger sizes. Thus, the size distribution 
at large sizes is very steep and g 4 is large. As the evolution proceeds, runaway growth 
produces more and more large objects from the ensemble of 100 km objects. With more 
objects at r > 500 km and fewer at 100 km, the power-law size distribution grows shallower 
and shallower (Fig. [2$; q<± decreases. In the late stages of oligarchic growth, planets reach 
their maximum radius. Because large objects mostly accrete mater ial from much smaller 



objects, the slope then converges on a characteristic value g 4 ~ 3 ( jGoldreich et al.l 12004 



Schlichting fc Saril 120111 ) . Although the growth of the largest objects is monotonic in time 
(Fig. H]), small number statistics produces the fluctuations around the steady decrease of g 4 
with time shown in the Figure. 

To examine the relationship between r max and g 4 , we consider calculations for each set 
of initial conditions (x m , r , g in j t , fi). For each set of calculations, we collect all pairs (r max , 
g 4 ) in all annuli at all times. We then bin these data in increments of 0.1 in log r max and 
calculate the median g 4 in each r max bin. This procedure yields the variation of median g 4 
as a function of r max and the initial conditions. This median provides a good representation 
of the typical slope: the inter-quartile range is roughly 2 when r max is small (~ 100 km) and 
roughly 0.5-1.0 when r max is large (> 1000 km). 

Fig. [6] shows results for the complete set of calculations at a = 15-75 AU with ro = 
1 km and the f w fragmentation parameters. For the ensemble of massive disks with x m = 3 
(violet points in the Figure), the median slope is large, g 4 ~ 10, when r max « 100-300 km. 
For r max ~ 200-2000 km, the median g 4 declines roughly monotonically with increasing r max . 
As r max grows from ~ 2000 km to ~ 5000 km, the median g 4 is roughly constant. 

The variation of g 4 with r max is sensitive to the initial disk mass. As the initial disk 
mass declines from x m = 3 to x m = 0.01, growth times are longer and longer. Due to 
the impact of gas drag and long-range stirring, small planetesimals in lower mass disks 
tend to have larger eccentricities at the onset of runaway growth than those in higher mass 
disks. Larger eccentricities reduce gravitational focusing factors, slowing growth and enabling 
stirring to grow faster. Because slower growth is more orderly, large objects have shallower 
size distributions. When r max m 100-300 km, lower mass disks thus have a smaller median 
g 4 than higher mass disks. Because growth is very slow for the smallest disk masses, all 
calculations with x m < 0.1 have roughly the same median g 4 « 4. 

For disks with x m > 0.3, gas drag (long-range stirring) is more (less) effective. With 
larger gravitational focusing factors in more massive disks, growth is relatively fast. Once a 
few 100 km objects form, a stronger runaway leads to a more pronounced decline of g 4 with 
r mnr . As from roughly 100 km to ~ 2000 km, the median g 4 declines roughly 

linearly with r max . Once r max > 2000 km, the median g 4 is roughly constant. In lower mass 
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disks, the median g 4 declines more slowly with increasing r max . Because low mass disks never 
produce objects with r max > 2000 km, the median g 4 is never constant with increasing r max . 

In addition to the clear dependence on x m , the median g 4 is also sensitive to the initial 
maximum planetesimal size r (Fig. [7]). For r max « 100-300 km and x m ~ 1-3, the median 
q4 steadily increases with larger ro, from g 4 8-9 for Tq = 1 km to g 4 ~ 10-11 for ro = 
10 km to g 4 ~ 13 for ro = 100 km. Independent of ro, the median g 4 declines steadily 
with larger r max . At the largest r max attained in our calculations, the median g 4 converges 
on an equilibrium slope ~ 3-4. Calculations with smaller r yield smaller values for this 
equilibrium slope. 

The behavior of the g 4 -r maa , relation is remarkably independent of semimajor axis (Table 
[2J Fig. [S]). In calculations at 30-150 AU, our 10 Gyr evolution times yield larger r max than 
the 1 Gyr evolution times for calculations at 15-75 AU. Otherwise, the median slope at small 
r max , the monotonic decline in slope at r max « 300-2000 km, and the roughly constant slope 
for r max > 2000 km closely follow the relations established at 15-75 AU. Our results show 
some indication for a smaller inter-quartile range in the median g 4 at larger semimajor axis. 
Confirming this trend requires a larger ensemble of calculations. 

Despite the insensitivity to semimajor axis, the median slope depends strongly on the 
initial size distribution (Tabled Fig. [9]). In calculations where the initial mass is distributed 
among planetesimals with sizes of 1 m to 1 km, the median g 4 is systematically smaller (~ 
1-2) throughout the evolution. This behavior is independent of r . All of these calculations 
begin with a substantial amount of mass in objects with sizes of 0.1 km or smaller. During 
runaway growth, collisional grinding removes some of this material from the grid. Compared 
to calculations where all of the initial mass is concentrated in the largest planetesimals, the 
largest objects in these calculations accrete material from a lower mass reservoir of small 
planetesimals which are easier to erode. Thus, the growth of the large objects is relatively 
slower and the median slope is smaller. Because calculations with larger ro have a smaller 
fraction of the initial mass in small planetesimals, this difference is smaller in calculations 
with larger r . 

To conclude this discussion, Fig. [TO] shows how the g 4 -r max . relation depends on the 
fragmentation parameters. For all x m , larger planets grow in disks composed of stronger 
planetesimals. Once the collisional cascade begins, collisions remove more mass from weak 
planetesimals than from strong planetesimals. At identical evolution times, disks composed 
of weak planetesimals therefore have less mass in bins with small planetesimals. With oli- 
garchs accreting most of their mass from small planetesimals, they do not grow as rapidly 
or as large when planetesimals are weak. 
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Despite the differences in r max , the median g 4 is fairly independent of the fragmentation 
parameters (Table[2]). In the most massive disks with x m = 0.33-3, calculations with stronger 
planetesimals have slightly larger g 4 than calculations with weak planetesimals: g 4 = 2.9-4.1 
for f w , g 4 = 3.2-4.3 for f s , and g 4 = 3.1-4.2 for f vs . These differences are smaller than 
the typical inter-quartile range of 0.5-1.0 when r max 1000 km. In lower mass disks with 
%m % 0.1, calculations with stronger planetesimals tend to have shallower slopes, but the 
inter-quartile ranges are still much larger than the differences. 

To summarize, Figs. IBHTU1 demonstrate a clear relation between the size of the largest 
object and the slope of the cumulative size distribution for r > 100 km. For a broad range of 
initial conditions and fragmentation parameters, g 4 is large (6-12) when r max w 300-500 km 
and small (2-5) when r max « 1000-3000 km. For any r max , disks with smaller masses have 
smaller g 4 than more massive disks. 

To understand how the the rest of the size distribution evolves with r max , we now 
consider the other three power-law exponents. We begin with the size distribution of 0.1- 
1 km objects in §3.21 and then discuss the evolution of 1-10 km (10-100 km) objects in §3.31 

3.2. The Size Distribution of 0.1-1 km Objects 

The slope of the cumulative size distribution for 0.1-1 km objects, qi, is a proxy for the 
evolution of objects smaller than vq. As the largest objects grow, they remove planetesimals 
from this size range. During runaway growth and the early stages of oligarchic growth, 
oligarchs accrete only a small fraction of the mass in 0.1-1 km objects (KB08, KB10). 
Thus, the slope of this size distribution (qi) should remain roughly constant with increasing 
r max - Once the collisional cascade begins, destructive collisions among small objects remove 
material from this range of sizes; destructive collisions among larger objects add material. 
Although the evolution of the slope then depends on the relative importance of destructive 
collisions among large and small objects, the collisional cascade eventually removes all of 
the mass in 0.1-1 km objects from every annulus (KB08, KB10). Thus, q\ should evolve 
dramatically during the cascade. 

Fig. HH shows the variation of the median q\ with r max for all calculations with the f w 
fragmentation parameters. Each panel illustrates results for a different value of r ; the lower 
(upper) set of points corresponds to calculations with initial n c oc r -0 5 (n c oc r~ 3 ). Colors 
indicate initial disk mass ranging from x m = 3 (violet) to x m = 0.01 (maroon) as in Fig. |6j 
The remarkable overlap of colored points demonstrates that the evolution of the median qi 
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is independent of initial disk mass. 

For calculations with r = 1 km, there are two main features of the evolution of gi with 
r max . As long as r max < 300-500 km, qi is constant. In the lower panel of Fig. [LTJ qi 
remains fixed at its initial value as the largest objects grow from sizes of a few km up to 
roughly 300 km. Thus, the size distribution of 0.1-1 km objects 'remembers' the initial size 
distribution throughout runaway and oligarchic growth. 

Once the collisional cascade begins, the general evolution of q\ is nearly independent of 
the initial size distribution and the initial size of the largest planetesimal. As the largest 
objects grow to sizes of ~ 1000 km, all size distributions begin to converge on q\ « 2.0-2.5 
(Table 3). After maintaining this limit for a small range in r max , qi approaches a limiting 
value of 0.0-0.5. As the cascade removes nearly all of the mass from an annulus, qi remains 
at this limiting value. 

The intermediate value of q\ ~ 2.0-2.5 results from collisional erosion (see the Appendix 
for a simple derivation). When the cascade begins, the ratio of the collision energy to the 
binding energy is Q c /Q*d ~ 0.1-1. In this regime, collisions produce a modest net reduction 
in the mass of the larger object. Together with the mass of the smaller object, this material 
is lost as debris and accumulates in much smaller mass bins. Smaller objects are easier 
to erode. Thus, the size distribution evolves from the growth limit of roughly 3-4 to an 
intermediate limit of roughly 2.0-2.5. 

The final limiting value of qi results from complete shattering (see the Appendix). As the 
cascade continues, the largest objects grow. Collision velocities among the smaller objects 
also grow. Once Q c /Q*d > 2-3, nearly all of the mass of two colliding objects goes into 
debris. Within the 0.1-1 km size range, all objects are susceptible t o shattering. When 



shatte ring of the small objects becomes important, qi approaches zero (IKenyon fc Bromley 



2004d). 



For calculations with r$ = 10-100 km, 0.1-1 km objects collide and merge into larger 
objects before the collisional cascade begins. Based on the evolution of the largest objects in 
Figs. I6T-TTU} growth tends to produce size distributions with q w 2-4. Thus, when the initial 
size distribution is steep (n c oc r -3 ), there is little evolution in q\ for r max < 300-500 km (Fig. 
[TTT upper sets of points in the middle and upper panels). When the initial size distribution 
is shallow (n c oc r~ a5 ), slow growth leads to a gradual evolution from g x « to q\ ~ 2-3 
(Fig. [TTJ lower sets of points in the middle and upper panels). 

Aside from this difference in the early evolution of qi, these calculations also approach 
the two limiting values earlier in their evolution (at smaller r max ). When ro = 10-100 km, 
there is initially less mass in the 0.1-1 km range; thus, the collisional cascade modifies qi 
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earlier. At the onset of the collisional cascade, calculations with — 1 km have most of 
their mass in 0.1-1 km planetesimals. Because it takes more collisional debris to change qi 
in these calculations, q\ begins to change when r max is larger. In these calculations, q\ also 
has a larger dispersion as a function of initial disk mass and r max . 

Analyzing calculations with other fragmentation parameters yields nearly identical re- 
sults (Table 3). The trend of the evolution is always independent of initial disk mass. For 
r = 1 km, qi always remains close to its initial value. When r = 10-100 km, qi remains 
roughly constant (when the initial q\ is larger than 2) or slowly evolves to q\ ~ 2-3 (when 
the initial q\ is smaller than 2). After the cascade commences, q\ first evolves to the canon- 
ical slope for erosive collisions of 2.0-2.5 (at r max « 1000 km) and then to the canonical 
shattering slope of roughly 0.0-0.5 (at r max « 3000-5000 km). Although it always occurs at 
r max ~ 500-1500 km, the transition from the initial q\ to the canonical q\ depends on the 
amount of mass initially in 0.1-1 km objects. Disks with less (more) mass in these objects 
make the transition at smaller (larger) r max and earlier (later) in evolution time. 



3.3. The Size Distribution of 1-10 km Objects 

Depending on the initial size of the largest object, the slope of the cumulative size 
distribution for 1-10 km objects, q2, serves two distinct roles. For calculations with r = 
1 km, the median q2 first measures changes in the slope as objects grow to sizes of 10 km 
and then tracks the changes in the size distribution of these objects throughout the rest 
of runaway /oligarchic growth and the collisional cascade. In these models, we expect the 
median q 2 first to decline with increasing r max (as in the evolution of q^ with r max ) and then 
to evolve to a slope characteristic of the collisional cascade. For calculations with larger 
objects, ro = 10-100 km, changes in the size distribution for 1-10 km objects should follow 
the changes in q\ observed for the 0.1-1 km objects: the median q 2 should evolve slowly to 
~ 2-3 during runaway /oligarchic growth and then evolve more rapidly to a smaller value 
during the collisional cascade. 

Fig. [12] shows the variation of the median q 2 with r max for all calculations with the f w 
fragmentation parameters. The format is identical to Fig. [TTJ with separate panels for each 
r , colors indicating the initial disk mass, and an upper and lower set of points corresponding 
to our two adopted initial size distributions. As in Fig. [TTJ the evolution of the median q 2 
as a function of r max is independent of the initial disk mass. 

In each panel, the evolution of q 2 follows expectation. In the lower panel (tq = 1 km), 
the median q 2 declines from ~ 5-6 when r max p=j 20 km to ~ 4 for r max > 100 km. This 
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evolution is independent of the initial size distribution. Once the collisional cascade begins, 
q 2 declines to ~ 2.0-2.5 as r max grows from 1000 km to 3000 km. At early times, the slope of 
the 1-10 km size distribution is similar to the slope of the large object size distribution (g 4 ); 
at late times, the slope evolves into one of the canonical slopes for the collisional cascade. 

For larger Tq, the evolution of q 2 tracks features of the evolution of q±. In the middle 
panel of Fig. [121 <?2 follows much of the qi evolution for 1 km objects (Fig. [TTJ lower panel). 
During runaway and oligarchic growth, q 2 remains constant with r max . Once the cascade 
begins, q 2 evolves to the canonical value for the early stages of the collisional cascade q 2 ~ 
2.0-2.5. As r max approaches 10 4 km, q 2 remains constant at roughly 2.5 instead of dropping 
to ~ 0. 

In the upper panel, the evolution of q 2 for r = 100 km similarly tracks the q\ evolution 
for r = 10 km. As the largest objects grow from ~ 100 km to ~ 1000 km, 1-10 km 
objects collide and merge into larger objects. Growth tends to evolve the slope to 2-4. For 
calculations with n c oc r -3 (upper set of points), q 2 does not change. When the initial 
size distribution is concentrated into large objects (n c oc r~°- 5 ), q 2 slowly approaches 2. 
Eventually, the collisional cascade begins; q 2 evolves to 2-2.5 and remains in this range as 
fmax reaches roughly 10 4 km. 

The ratio of the collision energy Q c to the binding energy Q* D allows us to understand 
the difference between the qi evolution for r = 1 km and the q 2 evolution for r = 10-100 km. 
When the collisional cascade begins, erosive collisions remove mass from the 0.1-1 km and the 
1-10 km size ranges. During this evolution, both qi and q 2 evolve to 2.0-2.5. Once 3000 km 
< r max < 10 4 km, collisions among smaller planetesimals produce shattering. However, 
the larger planetesimals remain in the erosive regime. As a result, objects are completely 
removed from the ensemble of 0.1-1 km objects (producing qi = 0.0-0.5) and slowly eroded 
among the ensemble of 1-10 km objects (leaving q 2 at 2.0-2.5). 

Results for calculations with the f s and f vs fragmentation parameters are nearly identical 
to those with the /„, fragmentation parameters. When collisions produce mergers instead of 
debris, q 2 evolves from its initial value to the canonical 'merger' limit of 2-4. When r is 
1 km, q 2 has no memory of the initial size distribution. Calculations with larger r remember 
the initial q 2 and maintain this slope throughout much of the evolution. Once the cascade 
begins, q 2 always evolves to ~ 2. Because 1-10 km objects are not completely shattered at 
any point in the collisional cascade, q 2 always remains larger than ~ 2 and never drops to 
0-1. 
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3.4. The Size Distribution of 10-100 km Objects 

Based on our results for q\ and q2, the evolution of q 3 , the slope for 10-100 km objects, 
should follow one of two paths. When the initial size of the largest planetesimal is 10 km 
or smaller, planetesimals must grow to fill the 10-100 km size bins. Growth produces a 
characteristic evolution of the slope from large values to ~ 2-4 (Fig. [6]). As the evolution 
makes the transition from growth to the collisional cascade, the slope evolves to a canonical 
value of roughly 2. Because 10-100 km objects are harder to shatter than 1-10 km objects, 
the slope produced by the cascade should never fall below 2. When planetesimals have 
initial sizes larger than 10 km, collisional growth requires a long time to erase the initial size 
distribution. Eventually, growth and the cascade modify the size distribution, which should 
approach the canonical value of 2. 

Fig. [13] confirms these expectations for calculations with the f w fragmentation param- 
eters. Following the format of Figs. [TlT[T2"| there is a panel for each r , colors indicating the 
range in x m , and sets of points for each initial size distribution. Once again, the evolution 
is nearly independent of the initial disk mass. 

For calculations with r = 1 km, growth initially produces large q 3 w 4-6. Shallower 
initial size distributions and more massive disks tend to yield larger q 3 . As objects grow, all 
calculations converge on a smaller range of q 3 , ~ 2.5-3.5, with smaller initial disk masses 
favoring smaller q 3 . For r max > 1000-2000 km, ^3 is independent of initial disk mass. 

When r = 10 km, growth leads to a smaller range of initial median slopes, ^3 ~ 5-6 for 
r max ~ 200 km. These median slopes are independent of the initial size distribution and the 
initial disk mass. As the largest objects grow to r max « 1000 km, q 3 reaches an approximate 
equilibrium at ~ 4. Once the cascade begins, this slope gradually declines to q% ~ 2. 

This evolution for (73 follows the evolution of q2 for calculations with ro = 1 km. In- 
dependent of the initial size distribution, growth produces a characteristic size distribution 
which slowly evolves to shallower slopes as the largest objects grow. Eventually, the cascade 
changes the size distribution and the median q 3 becomes smaller. In the middle panel of Fig. 
[T3"| this change occurs when r max > 3000 km. A similar, but less obvious, change in slope 
occurs for the tq = 1 km calculations in the lower panel. 

When r = 100 km, q 3 remembers the initial size distribution throughout most of the evo- 
lution (Fig. [131 upper panel). In these calculations, the initial size distribution is unaffected 
by runaway /oligarchic growth and the early stages of the collisional cascade. Throughout 
these phases of the evolution, growth is slow and collisions have little impact on large objects. 
Thus, the size distribution is unchanged. During the late stages of the cascade, occasional 
collisions among large objects produce some debris which fills all mass bins in the 10-100 km 
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size range. Independent of the initial size distribution, this debris drives the slope of the size 
distribution to the canonical value of q$ = 2.0-2.5. 

The evolution of q% is independent of the fragmentation parameters (Table 3). When the 
planetesimals are initially small, q% first evolves from large values (~ 4-6) to small values (~ 
3). Once the cascade begins, q% evolves to 2.0-2.5. In this size range, the collisional cascade 
has less impact than at smaller sizes. Thus, the transition from the q% ~ 3-4 characteristic of 
runaway/oligarchic growth to q% ~ 2.0-2.5 occurs at larger r max (and later evolution time). 



3.5. Summary of Size Distribution Evolution 

At 15-150 AU, the evolution of planetesimals follows a similar path at every distance 
from the central star. Initially, collisional damping and gas drag reduce the velocities of small 
objects. Among larger objects, low velocity collisions produce mergers instead of debris. 
Along with dynamical friction and viscous stirring, damping produces large gravitational 
focusing factors; this evolution leads to runaway growth of the largest planetesimals. As 
these protoplanets grow, they stir up the leftover smaller planetesimals. Eventually, collisions 
among the leftovers produce debris instead of mergers. A collisional cascade ensues, which 
grinds planetesimals with sizes of 1-3 km and smaller to dust. 

Despite the similarity in the evolutionary path, there is considerable diversity in the 
overall outcome. For identical starting conditions, the stochastic nature of coagulation leads 
to a broad range in the rate protoplanets grow, in the maximum radius of protoplanets, in 
the distribution of sizes as a function of time, and in the amount of mass lost through the 
collisional cascade. Although this diversity does not mask clear tre nds in the median outcome 
as a function of initial conditions (IKenyon &: Bromley! 120081 . |2010| ). the chaotic nature of the 
evolution prevents associating a single outcome with a single set of starting conditions. 

To quantify the overall evolution of the size distribution of icy planetesimals at 15- 
150 AU, we define four slopes, qi-q^, which measure the best-fitting slope of a power-law 
size distribution between 0.1-1 km (gi), 1-10 km (g 2 ), 10-100 km (g 3 ), and 100 km to r max 
(54). For all r max > 200 km, we derive the median slopes and the inter-quartile range as a 
function of r max . We use median qi-q^ and the inter-quartile range to track the evolution of 
the size distribution and to measure the range of outcomes as a function of initial conditions. 

Based on these slopes, we infer several general features of the evolution: 



• For all distances from the c entral star, the growth of icy planets leads to size distribu- 
tions with q « 2-4 (see also lKenyon fc Luulll999bl ; lKenyonll2002l ; iGoldreich et al.ll2004l ; 
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Schlichting fc Saril 120111 ) . This size distribution results from large objects accreting 
primarily small objects rather than other large objects. As largest objects grow, the 
slope evolves from g 4 pa 6-14 when r max pa 200 km to g 4 ~ 2-4 when r max > 1000 km. 

• The collisional cascade produces shallower slopes in the size distribution, q ^ 2.0-2.5. 
During the casc ade, several processes set the size distribution for objects with r > 
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mass from all objects and produces q ~ 2.0-2.5. Once the collisional cascade reaches 
its peak [r max > 3000 km), shattering of 0.1-1 km objects leads to q pa 0.0-1.0. 

The evolution of the large object size distribution is independent of radial distance 
from the central star and is generally independent of the fragmentation parameters. 
However, the evolution of g 4 is sensitive to the initial size of the largest planetesimal 
(ro) and the initial mass of the disk (x m ). Lower mass disks evolve over a smaller range 
of g 4 . Calculations with larger planetesimals tend to produce larger median g 4 . 

For any r max , the range in g 4 is fairly large. At small sizes [r max pa 200-500 km), the 
inter-quartile range about the median g 4 is roughly 2. This range drops to roughly 1 
at r max > 1000 km. For a typical g 4 pa 3, roughly 50% of calculations with the same 
initial conditions yield g 4 = 2.5-3.5. The rest have g 4 = 2-2.5 and g 4 = 3.5-5. 

Among the intermediate mass objects with r = 1-100 km, the evolution of the slope 
of the size distribution depends on the initial size of the largest planetesimal. 

When ro is small, growth produces objects with r = 1-100 km. Growth erases knowl- 
edge of the initial size distribution; the slope evolves to the standard 'growth slope,' 
q 2 pa q 3 pa 2-4. Once the collisional cascade begins, the slopes evolve to the canonical 
value of 2 for erosive collisions. 

When r is large, growth cannot change the initial slope of the size distribution. Thus, 
q 2 and q% remember the initial q until the onset of the collisional cascade. Once the 
cascade begins, erosive collisions lead to g 2 ~ 93 ~ 2.0-2.5. 
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4. APPLICATION TO THE TRANSNEPTUNIAN REGION 

To apply these results to the formation of TNOs, we place our calculations in the 
context of recent observations of TNOs and models for the formation of the solar system. 
The observations allow us to set r max and the slopes of the size distribution in various 
size ranges. Observations of young stars in nearby star-forming regions and theories of 
solar system formation constrain the initial properties of the protosolar disk and the likely 
epoch(s) of TNO formation. We begin with the theoretical issues and then summarize the 
direct observational constraints. 



4.1. Theoretical Constraints 



Current theories for the origin of the solar system broadly explain the nature of TNOs as 
a several step process. During the early evolution of the proto s olar disk, dust grains w ithin 
the disk grow into icy planetesim als fe.g.. IWeidenschil ling 20 06nChiang fc Youdinll2010) and 



then into i cy protoplanets (e.g., iLissauer 



1987 



2008 



Kenvon fc Luu 



2010) . Four protoplanets become gas giants (IPollack et al. 



999a: 



1996 



Kenvon fc Bromley 



Alibert et al. 



2005 



Ida fc Linll2008t ILissauer et al.ll2009t iBromley fc Kenyonll2011al ). Sometime after the gas gi- 
ants reach roughly their final masses, dynamical interactions among the gas giants, the much 
less massive icy protoplanets, and any leftover planetesimals result in the ej ection or migra- 
tion of icy objects into the transneptunian region and the Port Cloud (e.g.. iMalhotral 11993 



Hahn fc Malhotralll999l : Morbidelli et al.l2004l Tsiganis et al.l200a iGomes et al 



2005|). Once 



dynamical evolution begins, lower surface densities and larger orbi tal velocities among the 
TNOs reduce growth rates and promote collisio n al erosion (e.g ., 



2008 



Fraser 



200C 



Kenyon fc Luul Il999al: iKenvon fc B romlevj l2004cr lO'Brien et al.l l2005t Kenvon fc Bromley 



Stern fc Colwell 1997b 



9T Kenyon fc Bromley 2010 ). Thus, dynamical evolution ends the growth of 



TNOs. 

Despite the uncertainty in the theory, this picture constrains the initial mass of the 
protosolar disk and the likely range of formation times for TNOs. The protosolar disk 
probably has properties similar to the disks o bserved in nearby star- f orming regions such 
as Ophiuchus, Orion, and Taurus-Auriga (e.g.. IWilliams fc Ciezal 120111 ). Infrared and radio 
observat ions demonstrate that these systems lose their gas and dust on timescales of 1 - 



lOMyr flHaisch et al.ll200ll ; ICurrie et al.l 120091 : Kennedy & Kenvon! 120091 ; IWilliams & Cieza 
201 lh . For the solar system, these data constrain the initial disk mass, 0.001-0.1 M , and 
the timescale, 1-10 Myr, for the formation of gas giant planets. 



Dynamical calculations of the solar system place additional limits on the initial disk 
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mass and the formation time. Several results sugge st some TNOs formed at 20-25 AU and 
then migrated outwards with the giant planets (e.g.. lMalhotralll993[ Il995t lHahn k Malhotra 
19991 ; iLevison k Morbidellill2003l ; iMorbidelli et al.ll2008l ). To explain the orbital architecture 
of the giant planets, the TNOs, and the P ort cloud, these calculat i ons rely on a massive 
disk of solids with x m > 1 at 15-30 AU f Hahn k Malhotral Il999l : ILevison k Morbidelli 



2003; Morbidelli et al. 



20081 ). Although these calculations place little direct constraint on 



the formation time, Figs. [6T-TT01 demonstrate that more massive disks have larger g 4 at fixed 
r max- Thus, relying on a massive disk to enable outward migration places indirect constraints 
on any combination of formation time, r max , and 

Various lines of evidence associate migrating gas giant planets wi th an increased rate 



of cratering during the Late Heavy B ombardment (hereafter LHB; e.g.. ILevison et al.ll2001 



Gomes et al. 2005 ; Strom et al. 20051). Although co nstraints on the cratering history prior 
to this epoch are weak (e.g., IChapman et all 120071 ) . detailed analyses suggest the intense 
period of bombardment ended ~ 500-700 Myr after the format ion of the solar system (e.g., 
Tera et al.l 11974 ; lHartmann et al.1 12000| ; IStofner k Ryder! 120011 ) . It is unlikely that the gas 
giants migrated after the LHB. Thus, the end of the LHB sets an upper limit on the formation 
time for TNOs. 

For our analysis, we concentrate on the formation of TNOs before the gas giants dynam- 
ically rearrange the solar system. The initial conditions chosen for our calculations broadly 
match the range of initial disk masses and semimajor axes applicable to TNO formation. 
By following the evolution of TNOs for 1-10 Gyr, we cover the most likely range of forma- 
tion times, ~ 10-700 Myr. Thus, our calculations can provide clear tests of the coagulation 
picture for TNO formation. 



4.2. Observational Constraints 



Recent observations place excellent limits on the masses and radii of the largest TNOs. 

20 081) have at leas t one large object 



All of the major dynamical classes (e.g.. iGladman et al. 



with radius r fa 450-1300 km (e.g., ISheppard et al. 



2011 



Petit et al.l 120 111 ). All-sky surveys 



are now nearly complete to an R-band magnitude R fa 21, the approxima te brightness of 



an object with radius r fa 350 km and an albedo of 10% at 42 AU (e.g., ISchwamb et al. 
2010l ; ISheppard et al.l 1201 ll ). Thus, the sample of largest TNOs is probably complete at 40- 
50 AU. For large TNOs with binary companions, period m easurements yield reliable masses 
of roughly 1.5 x 10 24 g for Quaoar dFraser k Brownl 120101) to roug hly 1.5 x 10 25 g for Pluto 
( IBuie et al.ll2006l ) and 1.6 x 10 25 g for Eris (IBrown k Schallerll2007f ). In our calculations, this 
mass range corresponds to r fa 600-1400 km. Thus, a typical r max = 1000±400 km matches 
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the current population of the largest TNOs. 



For deep surveys with sufficient observations, the luminosity function - written as cr(m), 
the surface number density of objects brighter than magnitude m - is reasonably well-fit by 
a double power-law, 

1 + lO( a L- a s)(m eq -m ) 



a{m) 



a 



mo 



lQa L (m-m ) _|_ ^Q(a L -a s )(m eq -m )-a s (m-m ) 



(9) 



where a 



mo is a normalization factor, mo is a reference magnitude, m eq is the brightness where 

large sizes, and as is the 



slope of the size distribution at small 


sizes (e.g.. Petit et al. 


12006. 


2008 


20081: Fuentes et all 


20091; Fraser et al. 20081 Fraser & Kavelaars 


2009 


Schwamb et al. 2010; 


Sheppard et al. 


2011: Petit et al. 2011 


)■ 



uentes &: Holman 



Fraser et al.l 



2010 



Converting observations of the surface nu mber density to constraints on the size dis- 
tribution requires additional information (e.g., iFraser et al.l 120081 : iPetit et al.ll2008l ). If the 
cumulative size distribution and the distribution of albedos follow power laws with exponents 
q (for the size distribution) and (for the albedo), 



q = 5(« - 0/2) + 1 . 



(10) 



Altho ugh TNOs display a broad range of albedos (e.g.. IStansberry et al.ll2008l ; iBrucker et al. 



20091 ). current data suggest a roughly constant albedo of 2%-20% for r & 30-500 km an d 



a larger albedo of 20% to nearly 80% for r > 700 km (Fig. 3 of IStansberry et al.ll2008l ). 
Most objects comprising the luminosity function are faint; adopting « is therefore a 
reasonable assumption. The slope of the size distribution is then q = 5a + 1. 

Recent analyses indicate a range in the slope of the luminosity func tion at large sizes- 
Several deep surveys yield aj, ~ 0.70-0.76, implying q j 4.5-4.8 fo r J3 = (jFuentes fc Holman 



20081 : iFuentes et al.ll2009l : IFraser fc Kavelaarsll2009[ ). IFraser et al.l (120101 ) divided their sam- 
ple of the classical Kuiper belt into a set of cold objects with inclination i < 5° and a set of 
hot objects with i > 5°. Their results suggest aj, ~ 0.80 {qj, ~ 5) for the cold objects and 



aj, ~ 0.40 (qr, ~ 3) for the hot objects (see al s olLevison fc Ste 



11 



2001: Bernstein et al. 2004 



Fuentes &: Holmanll2008l : IFuentes et al.ll201ll ). ISheppard et al.l (120 111) also prefer a shallower 
slope for the large objects. 

Other observations also distinguish the hot and cold popu lations within the Kuiper 



belt. The largest cold KBOs are smal ler ( jLevison fc Sternll200ll ) and have a larger fraction 



of binary comp anions (iNoll et al.ll2008l) than the largest hot KBOs. The cold KBOs also have 



that t he cold KBOs are considerably redder than the hot KBOs (e.g. 



20001 : iTruiillo fc Brown! 120021 : IPeixinho et al.l 120041 : iGulbis et al.ll2006 



Tegler & Romanishin 


Peixinho et al. 


2008 
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This diversity suggests that the c old and hot KBOs formed in differen t locations or 



have had different collisional histories ( Levison &: Stern! l2001t iMalhotral Il993l ; iGomesI 12003 



Levison fc Morbidellill2003l ; iLevison et al.ll2008f ). In current dynamical models, the observed 



(a, e, i) distributions of the hot KBOs are a result of formation at 15-25 AU and subsequent 
migration and scattering during the migration of the giant pla nets. In this picture , the 
cold population formed at a larger distance, perhaps in situ (see iBatygin et al.l l2011aJ . and 
references therein). 

To assign r max for th e hot and cold populations, we use results from direct mass estim ates 



of individual KBOs (e.g.. lGrundv et al.ll2009l ; iFraser Sz Brownll2010l;lGrundv et al.ll201l[) and 



deep probes of the K BO luminosity function (e.g., IFraser fc Kavelaars 



2009 



Fraser et al. 



2010l ; iPetit et al.ll201ll ). These analyses consistently demonstrat e that the largest hot K BOs 



are systematically larger than the largest cold KBOs (see also ILevison fe Sternll200ll ). We 
infer a typical 'size ratio' > 1.5-2, suggesting r max w 1000 km for the hot objects and r max 
< 600 km for the cold objects. Although the true r max for the cold population might be 
less than 600 km, results for (g 4 , r max ) = (5, 600 km) are a reasonable proxy for any size in 
the 400-600 km (0.4 — 1.3 x 10 24 g) range. After deriving the constraints these observations 
place on the models, we will consider alternative models for a smaller r max = 100-200 km. 



4.3. Testing the Calculations 

To use the observational and theoretical constraints to test our calculations, we assume 
that our g 4 is equivalent to the observed q^. For our purposes, the g 4 derived for the 
cold KBOs is identical to the g 4 inferred from most deep surveys. Thus, we consider the 
implications of g 4 = 3 and g 4 = 5 for r max = 600-1400 km. After a basic application of 
our models to these data in §4.3. we derive a set of 'success fractions,' which measure the 
fraction of models that match the observations of g 4 and r max ( §4.3.2[) . This analysis yields 
two plausible models for the formation of the cold and hot TNOs: (i) ensembles of 1 km 
planetesimals in a massive disk with x m = 0.3-3 and (ii) ensembles of 10 km planetesimals in 
a less massive disk with x m = 0.03-0.3. To select the 'best' model, we consider the range of 
formation times at 20-50 AU ( §4.3.3[) . These results favor ensembles of 1 km planetesimals 
as the most likely starting point for the current populations of hot and cold TNOs. 
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4-3.1. Basic Tests 

The steep size distribution of the cold Kuiper belt objects (KBOs) suggests these objects 
form in a massive disk. Although many calculations yield g 4 ~ 5 when r max = 1000 km, the 
median g 4 for low mass disks is usually much smaller than 4 (Table |2]). Thus, more massive 
disks with x m ps 0.3-3 are strongly favored over low mass disks with x m < 0.3. Because 
the median g 4 also depends on the initial planetesimal size, the initial disk mass required to 
match the observations depends on Tq. Larger r$ favors lower mass disks. For all tq, models 
with x m ps 0.3-3 yield good matches to the observed g 4 . 

The shallow size distribution of the hot KBOs allows us to rule out a broader range 
of collision models. From Figs. I6HI01 and Table [21 calculations with 100 km planetesimals 
never achieve g 4 < 3.5 when r max = 1000 km. Thus, these models provide a poor match to 
the size distribution of large, high inclination KBOs. Results for calculations with 10 km 
planetesimals are a little more encouraging. Although calculations with a single planetesimal 
size fail, growth in a low mass disk (x m ps 0.03) with a range of planetesimal sizes yields g 4 
close to 3. 

For the hot KBOs, a broad range of calculations with 1 km planetesimals result in a good 
match to the observed g 4 with r max = 1000 km. Models with a range of initial planetesimal 
sizes have a smaller median g 4 than those with a single initial size; thus, these models provide 
a better match to the observations. Models with small initial disk masses achieve smaller 
q4 than models of massive disks. Current uncertainties in the slope allow a broad range of 
initial disk masses to match the observations. 



4-3.2. Success Fractions 

Having established rough limits on the initial disk mass and planetesimal size, we now 
quantify the range of calculations which match the observational constraints on r max and 
qA. To make these choices, we consider ensembles of results for each set of initial conditions 
and fragmentation parameters. From Table [fl there are three ensembles with r = 10 km 
and 100 km, and five ensembles with tq = 1 km. For each of these ensembles, we collect all 
pairs of (r max , g 4 ) and bin these data in increments of 0.1 in log r max . Within each r max bin, 
stochastic effects produce a broad range of g 4 . Thus, we identify the fraction of calculations 
within each bin that match an adopted target g 4 as a function of r max . To span the observed 
range, we consider two target slopes, g 4 = 3±0.25 and 5±0.25, for r max = 600 km, 1000 km, 
and 1400 km. 

To judge the uncertainties of the derived 'success fractions,' we examine the variation of 
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the success fraction as a function of r max . For a single set of initial conditions (x m , r , g^t, 
/i), the success fraction usually has an obvious trend (from larger to smaller success fraction 
or from smaller to larger success fraction) as a function of r max . We set the uncertainty in 
a single success fraction as the dispersion around this trend. Typically, this dispersion (and 
thus the adopted uncertainty) is 0.05-0.10. 

Fig. [TH shows results for a target = 5±0.25. The three sets of vertical panels plot 
the fraction of successful models for r max = 600 km (left panels), 1000 km (middle panels), 
and 1400 km (right panels). The three sets of horizontal panels plot the success fraction 
for calculations with ro = 1 km (lower panels), 10 km (middle panels), and 100 km (upper 
panels). Within each panel, colored lines show the variation of the success fraction as a 
function of the initial disk mass. Different colors indicate different initial mass distributions 
and fragmentation parameters as listed in the figure caption (see also Tabled]). 

Clearly, a broad range of calculations match the target slope at rates of roughly 20% 
to 50%. For r = 1 km and r max = 600 km (lower left panel), calculations with any set of 
initial conditions and low disk masses (x m < 0.1) fail to match the target slope. Models 
with larger initial disk masses (x m > 0.1) have larger success fractions. Within this range 
of x m , calculations with an initially large range of planetesimal sizes (green, orange, and 
magenta curves) are more successful than calculations with a single planetesimal size (violet 
and blue curves). As r max increases to 1000 km (lower middle panel) and to 1400 km (lower 
right panel), the success fraction systematically decreases. In these panels, calculations 
with a single planetesimal size are somewhat more successful than those with a range of 
planetesimal sizes. Once r max = 1400 km, however, few calculations successfully match the 
target slope. 

Calculations with larger r also match the target slope. For r = 10 km and 100 km, 
success fractions of 20% to 40% are common. Among the r = 10 km panels, there is a 
clear trend of larger success fraction with increasing initial disk mass. Although there is a 
small trend of more success for calculations with a single initial planetesimal size, the trend 
is much weaker than the trend with r max . Among the models with r = 100 km, less (more) 
massive disks yield larger success ratios when r max is small (large). When r max is small, 
calculations with a broad range of initial planetesimal sizes often match the data; results for 
a single initial planetesimal size rarely match the data. This trend reverses at large r max , 
when calculations with a single planetesimal size are often more successful. 

Fig. [15] repeats Fig. [TJ]for a target g 4 = 3±0.25. In these examples, calculations with 
1 km planetesimals are much more successful than those with large planetesimals. For any 
Tmax (any vertical column of panels in the figure), the success fractions decline monotonically 
from 20%-40% at r$ = 1 km to 0% for tq = 100 km. Although there is little trend in the 
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success fraction with r max for calculations with tq =10-100 km, there is a modest trend for 
r = 1 km: smaller r max requires disks with smaller masses to match the target slope. In all 
panels for r = 1 km, there is little trend with the initial size distribution of planetesimals 
or the fragmentation parameters: all of the initial conditions yield success ratios of 20% to 
40% for x m = 0.03-3. 

To refine these tests, we now include other constraints on the hot and cold populations. 
In the dynamical models, hot KBOs formed closer to the Sun and then migrated out to their 
current orbits. In the context of our coagulation models, objects closer to the Sun grow 
faster (Fig. H|) and produce shallower size distributions (Fig. [5]) over the same evolution 
time. Both of these features of collisional evolution agree with the observations: the hot 
KBOs are larger and have a shallower size distribution than the cold KBOs. 

To find a combination of starting conditions that match the observations, we consider 
models that match (g 4 , r max ) = (3, 1000 km) for the hot objects and (5, 600 km) for the cold 
objects. Using results in Figs. IT4T - [T5| we rule out any models starting with 100 km plan- 
etesimals. These calculations never match the shallow slope for the hot objects. However, 
two options with 1 km and 10 km planetesimals are equally likely: 

• 10 km planetesimals in a low mass disk (x m = 0.03-0.3) 

• 1 km planetesimals in a massive disk (x m = 0.3-3) 

In both cases, models with a single planetesimal size or a range of planetesimal sizes match 
the constraints reasonably well. Models with a single planetesimal match more often when 
r = 10 km; models with a range of sizes match more often when r = 1 km. 

4-3.3. The Formation Time 

To select the 'better' alternative among these two options, we add a final constraint 
- the formation time - to our analysis. The dynamical models require TNO formation 
at 20-40 AU sometime between 1-10 Myr and 500-700 Myr. The smaller of these two 
constraints is not useful. Only a few models produce TNOs in 10 Myr at 20 AU; none of 
our calculations produce TNOs in 10 Myr at 40-50 AU. The longer constraint is very useful. 
Many calculations require more than 1 Gyr to produce 600-1000 km objects at 40-50 AU. 
Thus, the dynamical models rule out these coagulation models. 

To select between the allowed combinations of r and x m , we consider the evolution time 
required to form 1000 km objects as a function of semimajor axis. Fig. [16] shows results 
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for calculations with the f w fragmentation parameters. Results for other calculations are 
similar. Individual panels plot results for r = 1 km (lower panels), r = 10 km (middle 
panels), and r = 100 km (upper panels) and two initial size distributions (left and right 
panels). Colors indicate the initial disk mass, ranging from x m = 3 (violet) to x m = 0.01 
(maroon) as in Fig. |6j 

These results follow standard expectations for coagulation calculatio ns. Formation 



times always scale with the initial radius of t he largest planetesimal (e.g.. iLissauerl Il987 



Kenyon fc Luulll998t iKenyon fc Bromley! 120101 ) . As r increases from 1 km to 100 km, the 



clusters of points shift to longer formation time s. In a disk with a surface density S oc 
i m qy 3 ^ 2 , the formation time scales as tf oc x~ 3 a 3 ( iGoldreich et al.l 120041 : Kenyon &: Bromley 



20081 ). As a increases, formation times in each panel grow rapidly. Similarly, models with 



smaller x m have longer formation times. 

The solid horizontal line in each panel indicates the approximate end of the LHB, 500- 
700 Myr. Applying the additional constraint that TNOs form at a ~ 20-50 AU yields 
constraints on all of our calculations: 

f 0.01 (r /l km) 1 / 2 (a/20 AU) 3 all planetesimal sizes, 
Xm ~ \ 0.03 (r /l km) 1 / 2 (a/20 AU) 3 one planetesimal size. ^ ^ 

When r = 100 km, these constraints are severe: in situ formation of TNOs at 40 AU requires 
x m > 1-3. Thus, calculations with very large planetesimals are viable only if the initial mass 
of the disk exceeds the minimum mass solar nebula. For smaller r , calculations with smaller 
x m yield 1000 km TNOs in 500 Myr at 40 AU. 

Requiring TNOs to form in 10 Myr places more stringent limits on the initial disk mass 
and the range of semimajor axes. In this case, TNOs can only form at a w 20 AU for massive 
disks (x m = 1-3) composed of 1-10 km planetesimals. Rapid formation times exclude all 
combinations of models with ro larger than 10 km. 

These conclusions are independent of the fragmentation parameters. When other initial 
conditions are held fixed, formation times for icy planets with r max fa 1000 km are similar 
for all three sets of fragmentation parameters. Although the size of the largest object varies 
with fi (Fig. [T0|) . the variation of with r max is similar for calculations with weak, strong, 
and very strong planetesimals. Thus, the formation time places strong constraints on the 
initial disk mass and the initial sizes of planetesimals. 

From eq. HU formation at 40-50 AU places the most severe constraints on the initial disk 
mass. For 1-10 km planetesimals, we adopt the less restrictive model with all planetesimal 
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sizes. Then: 

> f 0.1 — 0.2 1 km planetesimals, 
m ~ \ 0.3 — 0.6 10 km planetesimals. 

Because coagulation models with 10 km planetesimals match the observed range in g 4 and 
r max only for x m ~ 0.03-0.3, the combined constraints of the formation time, the size of the 
largest object, and the slope of the large object size distribution strongly favor calculations 
with 1 km planetesimals. 



4.4. Consequences of Small Cold KBOs 

Some observations suggest the largest cold K BOs have r mnT ~ 100-200 km instead of 



the r max = 400-600 km used in our analysis (e.g., iBrucker et al.ll2009l ; iGrundy et al.l 12009 



201ll ). To investigate how a smaller r max impacts our results, we now consider the range of 
coagulation models which yields an observed q = 5±0.25 for r max = 100-200 km. For this 
analysis, we use q$ as a proxy for the slope of the size distribution for cold KBOs smaller 
than r max . For r max = 100-200 km, g 4 provides a poor measure of the slope below 100 km; 
q% provides an accurate measure of the slope over 20-200 km. 

To identify coagulation models which match the observed q% for a smaller r max , we derive 
success fractions for all sets of initial conditions in Table [1] assuming a target q% = 5 ± 0.25 
and r max = 100-200 km. For calculations with r = 100 km and r max = 100-200 km, q 3 is 
identical to the initial slope (see Fig. [T3]). Because our adopted initial slope is g 3 = 0.17 or 
3.0, our calculations never yield the observed q 3 = 5. Thus, all of our calculations with r = 
100 km fail to match the target slope. 

Calculations with r$ = 100 km and a steeper initial size distribution probably cannot 
match the target slope. When the initial q$ > 3, most of the mass is in the smallest objects. 
For modest initial surface densities, x m > 0.1, growth of these small objects should produce 
shallower slopes with ^3 w 3-4. 

For all other initial conditions, we derive success fractions of 10% to 50%. When r max 
is much smaller than 500 km, the collisional cascade has no impact on the slopes of the 
size distributions. Thus, success fractions for a target (^3, r max ) = (5, 100-200 km) are 
independent of the fragmentation parameters. For calculations with all sizes or one size of 
planetesimals, there are weak trends of increasing success fractions in more massive disks. 
Lower mass disks produce slower growth, which tends to yield shallower size distributions. 
However, this trend is weak for small r max . 

Despite this lack of trends with f\ or x m , there are trends with the initial size and 
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size distribution of planetesimals. For r$ = 10 km, a larger fraction of our models match 
the target q 3 for r max = 150-200 km (40% to 50%) than for r max = 100-150 km (10% to 
30%). Calculations with r = 1 km yield smaller success fractions for r max = 150-200 km 
(10% to 35%) than for r max = 100-150 km (35% to 50%). In both cases, models with a 
single planetesimal size have somewhat smaller success fractions, but the differences are not 
significant. 

Thus, coagulation models match the observed q 3 = 5 for the cold KBOs when r max 
= 100-200 km. Calculations with r = 1 km (10 km) yield the target slope more often 
when r max = 100-150 km (150-200 km). With little preference for the disk mass or the 
fragmentation parameters, this analysis does not identify a clear preference for ro or x m . 



4.5. Summary of Results 

The most basic conclusion of our analysis is that the robust relation between r max and 
g 4 allows coagulation models to make clear predictions of (g 4 , r max ) for any deep survey of 
TNOs. Unambiguous relations between q 3 and r max also enable clear predictions for small 
TNOs. In principle, observations of (g 4 , r max ) or (g 3 , r max ) can falsify current coagulation 
models. In practice, current observations agree with the predictions. 

Using observational results from deep surveys of TNOs together with conclusions derived 
from dynamical models of the solar system, we place good constraints on coagulation models 
for the origin of TNOs throughout the solar system. Fig. [17] summarizes the main results of 
our analysis. As we add more observational and theoretical constraints to the analysis, the 
set of coagulation calculations which 'match' the constraints narrows. Thus, a complete set 
of observations for all TNO classes provides a rigorous test of coagulation models. 

To summarize current limits on the models, we step through the observational and 
theoretical constraints. For a single size distribution of TNOs, we derive constraints on 
coagulation models from g 4 and r max (top two rows of Fig. H7[) . 

• With a success rate of 20% to 50%, calculations with a broad range of initial conditions 
and fragmentation parameters can match a target combination (g 4 , r max ) which is 
consistent with observations of TNOs. 

• For a specific target (g 4 , r max ), the range of success rates allows us to rule out some 
ranges of initial conditions and fragmentation parameters. 

• When the target (g 4 , r max ) = (5, 1400 km), calculations with large planetesimals are 
more successful than those with small planetesimals (Fig. IT4"|) . 
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• Smaller target slopes, ~ 3, appear to rule out calculations with large planetesimals 
(Fig. H5]). 

• In all cases, models with large disk masses are more successful than those with small 
disk masses. Thus, models with massive disks composed of large planetesimals are 
much more likely to match the target than low mass disks composed of small planetes- 
imals. 

Adding in the dynamical constraint that the hot and cold populations originate in 
different locations of the protosolar disk yields more severe limits on the models. Requiring 
(<?4, r max ) = (3, 1000 km) for the hot KBOs and (q^, r max ) = (5, 400-600 km) for the cold 
KBOs (third row of Fig. H7|) leads to these conclusions. 

• Models with 100 km planetesimals cannot match the data for any set of initial condi- 
tions. 

• Models with 10 km planetesimals match the data for low mass disks, x m = 0.03-0.3. 

• Models with 1 km planetesimals match the data for massve disks, x m = 0.3-3. 

Finally, using the constraint on formation time eliminates models with 10 km planetes- 
imals. Forming TNOs in 500-700 Myr at 20-50 AU requires massive disks with x m = 0.3-3 
(fourth row of Fig. fTT|) . Thus, we strongly favor models starting with 1 km planetesimals. 

Adopting a different observational constraint, r max = 100-200 km, for the cold KBOs 
results in similar conclusions. The observed slope for the hot population, g 4 = 3, drives us 
to models with small planetesimals in massive disks. 



4.6. Discussion 

Our analysis is the first attempt to apply multiannulus coagulation calculations to mod- 
ern observations of TNOs. The results indicate that TNOs form within a massive disk com- 
posed of small, 1 km planetesimals. This conclusion is similar to recent results on the origin 
of asteroids, where an initial population of 0.1 km plan etesimals seems suffic ient to explain 



the size distribution of asteroids for sizes of 10-100 km (jWeidenschillingl 1201 ll . for a different 
conclusion, see Morbidelli et al. 2009). Although we did not consider formation of TNOs 
from 0.1 km planetesimals, calculations with smaller planetesimals should yield similar re- 
sults for the size distributions at 10-1000 km. Thus, current data are also consistent with 
TNO formation from 0.1 km planetesimals. 
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The similarity of the size distributio ns for the TNOs and the tro jans of Jupiter (IJewitt et al 



2000l ; IJedicke et al.ll2002l ) and Neptune (jSheppard fc Truj illol |2 1 Oh suggest a common origin 



in a massive disk composed of 1-10 km planetesim als (for a discussion of the dynamical ori- 
gin of these populations, see iMorbidelli et al.ll2009l ). Unless the initial surface density of the 
disk is small (x m < 0.03-0.1; see eq. [IT]), collisional evolution erases the signatures of many 
initial conditions on 1 0-500 Myr timescales. Dynamical models for the orbital architecture 
require massive disks ( iMorbidelli et al.l 120081 ) . During the dynamical re- arrangement of the 
giant planets, the large relative velocities among planetesimals leads to a collisional cascade 
which modifies the size distribution at a w 15-30 AU considerably for r < 100 km (Figs. 

Q3HI3D- 

The main alternative to this picture is that the observed size distributions for the 
cold KBOs and the trojans of Jupiter and Neptune are 'primordial,' relatively unchanged 
since the initial process th at formed planetesimals from mm-sized or smaller particles (e.g., 
Sheppard fc Truj illol |2 1 01 ) . To evaluate this possibility, we isolate coagulation calculations 
where an initial size distribution of 100 km and smaller objects produces 110 km or smaller 
objects prior to the end of LHB. For models with a single size or a range of sizes, calculations 
with x m < 0.0003 (a/40 AU) 3 satisfy this constraint. Models with larger x m yield objects 
larger than 110 km. For primordial objects formed in situ at a « 20-30 AU, this limit 
implies an i nitial mass of 0.003 M m , roughly three times larger than the current mass in cold 
KBOs (e.g.. lBrucker et al.ll2009l ). Although coagulation models do not directly preclude that 



cold KBOs and trojans form in a very low mass region of the disk, it seems unlikely that 
dynamical mechanisms can popul ate three distinct volumes of the solar system with so little 
reduct i on in the total m ass (e.g., IMorbidelli et al.l 120051 ; iLevison et al.ll2008l ; iBatygin et al. 
201 lbl : IWolff et al.l 1201 if ). Thus, we conclude that the observed size distributions are not 
primordial. 

Comparisons between our calculations and the current size distributions of TNOs make 
several important assumptions. We assume that the hot and cold populations in the Kuiper 
belt are representative of all TNOs. Curr ent data do not contradict this assertion (e.g., 
Sheppard &: Truj illol 1 2 1 0t iPetit et al.ll201ll ). Data from ongoing and planned all-sky surveys 
will eventually yield robust estimates for (g 4 , r max ) for all dynamical classes among the TNOs 
and will likely improve estimates for the trojans of Jupiter/Neptune. If other dynamical 
classes of TNOs or the trojans of Jupiter/Neptune have different combinations of (54, r max ) 
than analyzed here, our approach provides a framework to test coagulation models with 
these new data. 

Detailed tests of our coagulation calculations also assume that large-scale rearrangement 
of the planets in the solar system has little impact on relations between r max and g 4 . We do 
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not expect this evolution to produce significant changes in q±. Before dynamical interactions 
with the gas giants become important, gravitational focusing factors are already small. Col- 
lision rates are in the orderly regime, where all objects collide at comparable rates. Thus, q^ 
should not evolve considerably once large-scale dynamical evolution is important. 

Throughout large-scale dynamical evolution, we expect significant changes to q\-qz as a 
result of collisional evolution. Because collisions remain in the orderly regime, these changes 
should follow the paths outlined in our discussion to Figs. [TTHTBl Depending on the timing 
of the start of the dynamical evolution, q 2 and q 3 for dynamically hot populations should 
evolve to roughly 2 or remain roughly constant at 2. If collision rates are large enough, q\ 
should evolve from roughly 2 to roughly 0. 

Some change in r max during the epoch of large-scale dynamical evolution seems likely. 
Because collisions occur at high velocities, it is unlikely that r max increases. Low collision 
rates im ply little collisional e rosion. However, the occasional high velocity, 'hit-and-run' 
collision ( lAsphaug et al.l 120061 ) could remove the outer shell of material from a large object 
and leave the denser cor e unchanged. This me chanism is one way to explain the relatively 
high density of Quaoar (IFraser fc Brown! |2010| ). Dynamical evolution can also eject a few 
of the largest TNOs from the solar system. Because the number of large objects is always 
small, a few ejections can reduce r max and lead to variations in r max among different TNO 
populations. 

Deriving the outcome of the dynamical evolution requires codes which combine statis- 
tical calculations of sm all objects as in this paper w ith the n-body calculations of gas giant 
planet formation as in iBromley &: Kenvonl (1201 laf) and calcu l ations of the migration an d 
scat tering of small obje c ts as i n ICharnoz fc Morbidellil (120071 ) , iLevison et al.l (120081 . 120101 ) , 
and lBromley fc Kenvonl (j2011bl ). Following the relative timing of the formation of gas giants 
and TNOs within the same disk will provide better limits for r max , g 4 , and possible locations 
for the formation of TNOs within the protosolar disk. Combining collisional growth with 
scattering will yield better predictions for the break radius and the slope of the size distri- 
bution below the break. Our encouraging results from pure coagulation models suggest a 
more complete model will enable a better understanding for the origin of TNOs and other 
large objects throughout the solar system. 



5. SUMMARY 



To develop a framework to test coagulation models for the formation of TNOs, we 
have analyzed a set of calculations in grids extending from 15-75 AU and from 30-150 AU. 
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The calculations consider a broad range of initial conditions for (i) the disk mass (x m = 
0.01-3 times the mass of the minimum mass solar nebula), (ii) the initial size of the largest 
planetesimals (r = 1, 10, and 100 km), (iii) the initial size distribution of planetesimals 
(n c oc r~ qinit , with q init = 0.5 and 3), and the fragmentation parameters (weak, strong, and 
very strong). Table [1] summarizes the number of calculations with each set of parameters. 

To facilitate comparisons between the calculations and observations, we focus on several 
observable features of the models. The size of the largest object in a region, r max , is easily 
derived in the models; in any survey, the brightest TNOs are usually the largest and most 
massive. If the albedos of TNOs are well-known, the slope of the size distribution follows 
from observations of the TNO luminosity function. Here, we define four quantities, qi-q^, 
which measure the slope of the size distribution over a decade in radius. Tables [2H5] list 
median values for these slopes as a function of r max for each set of initial conditions. 

The evolution of q\-q± as a function of r max is equivalent to the evolution as a function 
of time. Thus, we treat the evolution of icy planets at all distances from the Sun on the same 
footing and derive more robust relations between the observables and the initial conditions. 
For the theoretical calculations in this paper, we infer the following conclusions (see also 

MM-- 



For identical initial conditions, small variations in collision probabilities - due to dif- 
ferent random number seeds - produce a broad range of outcomes for g 4 , r max , and 
other observables. 



• As in iKenyon fc Bromley! (120081 . 120101 ) . the formation time for objects with r = r max 
> 300-1000 km depends inversely on the initial disk mass and with t he cube of the 
distance from the central star, tt oc x~ 3 a 3 (see also Safronov 1969; Lissauer 1987: 



Goldreich et al.l l2004h . 



• Coagulation models yield a robust correlation between r max and g 4 as a function of the 
initial conditions in the protosolar disk. This correlation allows unambiguous tests of 
the models. 



• When 1000 km protoplanets grow from smaller objects, the evolution naturally leads 
to a roughly power-law size distribution for objects with r > 100 km. The typical 
slope for the cumulative s i ze distribution is 04 ~ 2 -4 (see also IKenyon fc Luul Il999b 



Kenyon fc Bromley! l2004d ; ISchlichting fc Sarill201ll ). This slope is independent of the 



distance from the Sun and the initial size distribution. Although the slope depends 
weakly on the fragmentation parameters, it is very sensitive to x m and r$. Larger x m 
and larger r$ yield steeper size distributions with larger g 4 . 
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• For planetesimal sizes smaller than the initial size of the largest planetesimal (ro), the 
size distribution retains some memory of the initial size distribution. However, this 
memory is limited to timescales before the onset of the collisional cascade. Once the 
cascade begins, the size distribution depends on the evolution of the largest objects 
and the fragmentation parameters for the smallest objects. 

• Once protoplanets reach radii of 500-1500 km, high velocity collisions among leftover 
small planetesimals start a collisional cascade which grinds many of the leftovers to 
dust. As the cascade develops, the slopes of the size distribution for 0.1-100 km 
objects evolve to a standard value, q ~ 2.0-2.5. Eventually, the largest protoplanets 
reach sizes of 2000 km or larger. High velocity collisions among the leftovers then 
completely shatter 1 km and smaller objects. The slope of the size distribution for 
0.1-1 km objects then evolves to much smaller values, q\ ~ 0. Slopes for larger size 
ranges stay roughly constant at q% ~ q% ~ 2. 



To test the ability of coagulation models to match the observations, we focus our analysis 
on the relation between g 4 and r max for the hot and cold populations of KBOs. The q^ 
for the cold population and the r max for the hot population are representative of the entire 
ensemble of TNOs. Treating the two populations separately allows us to test the models as 
much as possible with current data. 

Our tests are very encouraging. A broad range of coagulation calculations match the 
typical slope, q 4 = 5, for large objects when r max = 400-1400 km (4 x 10 23 — 1.5 x 10 25 g). 
Large slopes generally favor calculations starting with 10-100 km planetesimals. However, 
smaller r max favors smaller planetesimals. The observed slope of the hot objects, q 4 = 3, 
completely rules out models starting with 100 km planetesimals. If TNOs must form prior 
to the LHB at 500-700 Myr, formation models for the hot objects require ensembles of 1 km 
planetesimals in a massive disk (x m = 0.3-3). If the largest cold objects have r « 400- 
600 km (4 x 10 23 — 1.5 x 10 24 g), the LHB constraint similarly requires calculations with 
1 km planetesimals. 

Combining the observations of the hot and cold populations with the constraints that 
they formed in different locations of the protosolar disk prior to 500-700 Myr leaves us with 
one set of starting conditions, a massive disk with 1 km planetesimals. This conclusion agrees 
with our previous result that t he basic building blocks of debris disks are 1 km planetesi- 



mals (IKenyon fc Bromley! |2010| ). Alt hough these analyses d o not consider the possibility of 



0.1 km or smaller planetesimals as in IWeidenschillingi ( 120111 ). the results clearly favor small 



planetesimals over larger planetesimals with radii of 10-100 km. 



Currently, it is difficult to make robust tests of predictions for Many surveys sug- 
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gest the TNO size distribution is a double power law with a break at R fa 25-26 (IFraser et al. 
20081 : IFraser fc Kavelaarsl bood : iFuentes et aDl2009h . If TNOs have a typical albedo of 6%, 
the break in the size distribution is at r fa 30 km. Although the uncertainties in the ob- 
servations below this break are large, the slope at small sizes is most likely shallower than 
at large sizes, with qs ~ 2 a reasonable estimate. The coagulation calculations described in 
this paper can match one of these two c onstraints. These calculat i ons produce a break in 



the size distribution at 1-3 km (see also iKenyon fc Bromley! l2004d ; iGil-Hutton et al.l 12009 



Fraserll2009l ). Below this break, the collisional cascade always produces a slope of roughly 2. 
Thus, our models in this paper match the slope below the break but not the location of the 
break. 



Kenyon &: Bromley! ( I2004cf ) demonstrate that stirring from th e largest nearby planet 
sets the 'break radius' (see also IGil-Hutton et al.ll2009l ; IFraser! 120091 ) . For TNOs, continued 
stirring by Neptune at 30 AU yields a break at 10-30 km, close to the lower limit inferred 
from current observations. Given the uncertainties in the luminosity function at R = 25-27 
and the albedos for such TNO s, coagulation models provide a satisfactory match to the data 
(see also iNesvorny et al.l 1201 ll) . 



Improving these tests requires additional data. Among the more massive TNOs, deriving 
q4 and r max for all of the dynamical classes will place better limits on the ability of coagulation 
and dynamical calculations to provide a unique model for the origin of the giant planets 
and TNOs. As deep surveys probe larger volumes of transneptunian space, results for the 
observational analogs to qi-q^ will test the ability of coagulation models to explain size 
distributions for each class of TNOs. For the entire ensemble of TNOs, coagulation models 
probably have fewer free parameters than the set of observational constraints. Thus, TNOs 
should provide robust tests of coagulation models for planet formation. 
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A. APPENDIX 

As oligarchic growth proceeds, cratering collisions gradually dominate physical inter- 
actions among 0.1-10 km planetesimals. In this regime, collisions between equal mass ob- 
jects produce the most debris. For a power-law cumulative size distribution with q > 3, 
small planetesimals are more likely to collide wi th one another than with large oligarchs 



( iKenyon fc Bromley! l2004at iGoldreich et al.ll2004j). Thus, cratering leads to a g radual ero 



sion of small planetesimals (e.g., iTanaka fc Idalll996l ; iKobayashi fc Tanakall2010l ). 



During this erosive phase, collisions between equal mass objects produce a single object 
with mass comparable to the target and copious debris. For these objects, the slope of the 
size distribution then depends on the rate collisions convert planetesimals into debris and the 
rate of debris production over a range of sizes. To estimate the rate planetesimals are lost, 
we consider the collision rate, navf g . During oligarchic growth, small planetesimals have 
similar orbital eccentricities and f g ~ 1. For collisions between equal mass objects, the rate 
small planetesimals are destroyed is dN/dt oc N 2 r 2 , where N is the number of planetesimals 
in a size range and r is a typical radius. Setting N oc r~ q , 

d log N/dt w r 2 ~ q . (Al) 

For an ensemble of planetesimals, small objects are destroyed relatively more (less) often 
when q > 2 (q < 2). More destructive collisions among small objects reduces the slope; 
fewer destructive collisions raises the slope. Without any debris, this process leads to an 
equilibrium size distribution with q w 2. In our simulations, debris has a power-law cu- 
mulative size distribution with qdebns — 2.5. Thus, cratering produces a power-law size 
distribution with q cr 2.0-2.5, as observed in our calculations. 

As the largest objects continue to grow, collisions among small objects become more and 
more catastrophic. For the fragmentation parameters we adopt, Q* D monotonically increases 
with r for > 0.1 km. In this range, small objects fragment catastrophically before large 
objects. Thus, catastrophic fragmentation reduces the equilibrium slope from the cratering 
slope of 2.0-2.5. 

To provide a rough estimate of the equilibrium slope for catastrophic collisions, we 
consider a collision rate where the fraction of completely destructive collisions is a decreasing 
function of r, e.g., fd oc r~ p . The loss rate of planetesimals is then, 

d log N/dt « f d r 2 - q ^r 2 ~ p - q . (A2) 

From the expression for Q* D (eq. [6]), p ~ 1.5. In the absense of debris production, catas- 
trophic collisions lead to a power-law size distribution with q c d ~ 0.5, close to the q ~ 0-1 
derived in our simulations. 
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Table 1. Grid of Calculations a with X oc a 3 / 2 



Initial size of largest object (ro) 
x m 1 km 10 km 100 km 1 km 10 km 100 km 1 km 10 km 100 km 1 km 1 km 



0.01 


7 


7 


7 


7 


14 


7 


7 


8 


10 


14 


15 


0.03 


7 


7 


7 


7 


7 


7 


7 


8 


9 


16 


15 


0.10 


7 


7 


7 


7 


12 


7 


7 


7 


8 


14 


15 


0.33 


7 


7 


7 


8 


14 


9 


7 


7 


15 


19 


17 


1.00 


7 


7 


7 


8 


16 


7 


14 


7 


13 


15 


16 


3.00 


7 


7 


7 


7 


7 


7 


13 


7 


13 


18 


15 


Notes 


b, e, g 


b, e, g 


b, e, g 


c, e, g 


c, e, g 


c, e, g 


c, f, g 


c, f, g 


c, f, g 


c, f, h 


d, f, 



a Number of independent calculations for each combination of x m and ro 
b a = 15-75 AU, t max = 1 Gyr, new calculations for this paper 
c a = 30-150 AU, t rnax 
d a = 30-150 AU, t max 

G n rv nr> 0-5 



10 Gyr, calculations from lKenvon fe Bromlevl (|2008l . 2010) 
10 Gyr, new calculations for this paper 



n c oc r 



l f w fragmentation parameters ([Leinhardt fc Stewartl l2009) 



h f s fragmentation parameters (jBenz fc Asphaudll999l ) 



1 f vs fragmentation parameters (|Davis et al.l ll985) 
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Table 2. Median values for q 4 when r ma:c = 1000 km 



Initial Disk Mass (x m ) 
Parameter 0.01 0.03 0.10 0.33 1.00 3.00 Notes 



94,3 1-78 2.08 3.62 4.13 5.06 5.30 a = 15-75 AU, r = 1 km, n c oc r" 5 , f w 

94,3 ••• 3.86 4.39 4.69 5.61 5.60 a = 15-75 AU, r = 10 km, n c oc r -0 ' 5 , /w 

94,3 5.25 6.27 6.53 a = 15-75 AU, r = 100 km, n c oc r" - 5 , jV 

94,3 1-83 2.79 3.64 5.48 6.31 6.22 a = 30-150 AU, r = 1 km, n c oc r -0 ' 5 , fw 

94,3 ••• 3.72 4.50 4.94 5.65 5.53 a = 30-150 AU, r = 10 km, n c oc r -0 ' 5 , 

94,3 5.03 5.42 6.61 6.93 a = 30-150 AU, r = 100 km, n c oc r" 5 , / w 

94,3 ••• 2.74 2.94 2.93 3.53 4.08 a = 30-150 AU, r a = 1 km, n c oc r~ 3 , / w 

94,3 • • • 2.72 3.73 4.68 5.18 5.70 a = 30-150 AU, r = 10 km, n c oc r~ 3 , / w 

94,3 ••• 3.52 4.03 4.74 5.12 5.63 a = 30-150 AU, r = 100 km, n c oc r~ 3 , /w 

94,3 1-85 2.18 2.58 3.17 3.55 4.25 a = 30-150 AU, r = 1 km, n c oc r" 3 , f s 

94,3 ••• 2.52 2.98 3.08 3.81 4.24 a = 30-150 AU, r = 1 km, n c oc r" 3 , /„ 5 



Table 3. Median values for q± when r max = 1000 km 



Initial Disk Mass (x m ) 



Parameter 


0.01 


0.03 


0.10 


0.33 


1.00 


3.00 








Notes 


91,3 


1.97 


2.00 


1.28 


1.13 


1.06 


1.29 


a 




15-75 AU, r = 


= 1 km, n c oc r~ ' 5 , f w 


91,3 




2.23 


2.23 


2.24 


2.26 


2.26 


a 




15-75 AU, r = 


- 10 km, n c oc r -0,5 , fw 


91,3 








1.99 


2.06 


1.91 


a 




15-75 AU, r = 


= 100 km, n c oc r -0 - 5 , f w 


91,3 


2.04 


1.56 


1.34 


1.19 


1.25 


1.40 


a 




30-150 AU, r 


= 1 km, n c oc r -0 - 5 , fw 


91,3 




2.25 


2.25 


2.25 


2.25 


2.26 


a 




30-150 AU, r 


= 10 km, n c oc r -0 ' 5 , fw 


91,3 






2.06 


1.98 


2.16 


1.78 


a 




30-150 AU, r 


= 100 km, n c oc r -0 - 5 , /iy 


91,3 




2.42 


2.17 


2.15 


2.20 


2.34 


a 




30-150 AU, r 


= 1 km, n c oc r -3 , /^y 


91,3 




2.18 


2.28 


2.39 


2.18 


2.20 


a 




30-150 AU, r 


= 10 km, n c oc r~ 3 , /vi/ 


91,3 




2.14 


2.14 


2.12 


2.11 


2.10 


a 




30-150 AU, r 


= 100 km, n c oc r~ 3 , /vi/ 


91,3 


2.03 


2.26 


2.22 


2.41 


2.43 


2.41 


a 




30-150 AU, r 


= 1 km, n c oc r" 3 , f s 


91,3 




2.39 


2.51 


2.29 


2.31 


2.37 


a 




30-150 AU, r 


= 1 km, n c oc r -3 , /„ a 
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Table 4. Median values for q 2 when r max 


= 1000 km 


Initial Disk Mass (x m ) 




Parameter 0.01 0.03 0.10 0.33 1.00 3.00 


Notes 



92,3 


3.25 


3.18 


3.63 


3.64 


3.61 


3.51 


a 




15-75 AU, r = 


= 1 km, n c oc r~ ' 5 , f w 


92,3 




0.94 


0.89 


0.89 


0.91 


1.02 


a 




15-75 AU, r = 


= 10 km, n c oc r -0,5 , fw 


92,3 








2.27 


2.24 


2.24 


a 




15-75 AU, r = 


= 100 km, n c oc r~ - 5 , fw 


92,3 


3.30 


3.59 


3.61 


3.66 


3.50 


3.35 


a 




30-150 AU, r 


= 1 km, n c oc r -0 - 5 , fw 


92,3 




0.95 


0.90 


0.93 


1.13 


1.05 


a 




30-150 AU, r 


= 10 km, n c oc r~ - 5 , fw 


92,3 






2.25 


2.26 


2.22 


2.25 


a 




30-150 AU, r 


= 100 km, n c oc r -0 5 , f w 


92,3 




4.03 


3.70 


3.72 


4.00 


3.83 


a 




30-150 AU, r 


= 1 km, n c oc r -3 , f w 


92,3 




2.96 


3.12 


3.14 


3.14 


3.15 


a 




30-150 AU, r 


= 10 km, n c oc r -3 , fw 


92,3 




2.56 


2.61 


2.69 


2.70 


2.73 


a 




30-150 AU, r 


= 100 km, n c oc r~ 3 , 


92,3 


4.44 


4.11 


4.00 


3.94 


4.00 


3.92 


a 




30-150 AU, r 


= 1 km, n c oc r~ 3 , / s 


92,3 




4.15 


4.00 


3.97 


4.05 


4.08 


a 




30-150 AU, r 


= 1 km, n c oc r -3 , /„ a 



Table 5. Median values for q 3 when r max = 1000 km 



Initial Disk Mass (x m ) 



Parameter 


0.01 


0.03 


0.10 


0.33 


1.00 


3.00 








Notes 


93,3 


4.17 


4.14 


3.85 


3.82 


3.71 


3.70 


a 




15-75 AU, r = 


= 1 km, n c oc r~ ' 5 , f w 


93,3 




4.11 


4.05 


4.06 


4.05 


3.98 


a 




15-75 AU, r = 


- 10 km, n c oc r -0,5 , fw 


93,3 








0.97 


0.93 


0.95 


a 




15-75 AU, r = 


= 100 km, n c oc r -0 - 5 , f w 


93,3 


3.87 


3.72 


3.67 


3.64 


3.70 


3.69 


a 




30-150 AU, r 


= 1 km, n c oc r -0 - 5 , fw 


93,3 




4.03 


4.06 


4.02 


3.96 


3.94 


a 




30-150 AU, r 


= 10 km, n c oc r -0 - 5 , fw 


93,3 






0.93 


0.97 


0.90 


0.92 


a 




30-150 AU, r 


= 100 km, n c oc r~°' 5 , / w 


93,3 




3.18 


3.43 


3.19 


2.93 


3.40 


a 




30-150 AU, r 


= 1 km, n c oc r -3 , /^y 


93,3 




3.87 


3.86 


3.82 


3.67 


3.71 


a 




30-150 AU, r 


= 10 km, n c oc r -3 , fw 


93,3 




3.16 


3.15 


3.16 


3.16 


3.17 


a 




30-150 AU, r 


= 100 km, n c oc r~ 3 , /vi/ 


93,3 


2.91 


3.20 


3.34 


3.36 


3.25 


3.35 


a 




30-150 AU, r 


= 1 km, n c oc r" 3 , / s 


93,3 




3.44 


3.84 


3.69 


3.59 


3.45 


a 




30-150 AU, r 


= 1 km, n c oc r~ 3 , /„ a 
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Fig. 1. — Variation of Q* D with radius for the three sets of fragmentation parameters adopted 
in the text. For the strong and very strong planetesimals, the solid (Q& = 10), dashed 
(Qb = 10 3 ) and dot-dashed curves (Qb = 10 5 ) indicate a range of parameters for the strength 
regime. 
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Fig. 2. — Time evolution of the cumulative number distribution for 8 annuli (a = 27-32 AU) 
in a disk with initial surface density distribution £ = 0.18 x m (a/30 AU)~ 3//2 g cm -2 , x m = 
1, and the f w fragmentation parameters. The ensemble of planetesimals has an initial cumu- 
lative size distribution n c oc r -0 5 with a maximum radius of 1 km; thus, most of the initial 
mass is in the largest objects. During runaway growth at 1-10 Myr, the largest protoplanets 
grow from ~ 10 km to ~ 1000 km. At the same time, the slope of the size distribution 
becomes shallower and shallower. At late times (t > 100 Myr), the size distribution is a 
fairly smooth power law from ~ 3-10 km to ~ 3000 km. 
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Fig. 3. — Definitions of the four slope parameters efined in the text. The black curve repeats 
the 300 Myr size distribution from Fig. [2j The slopes qi-q^are the slopes of the four line 
segments in the figure. 
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log Time (yr) 



Fig. 4. — Time evolution of the size of the largest object in 4 sets of eight annuli for the 
calculation described in Fig. [2J After a short period of slow growth, objects rapidly grow 
from ~ 10 km to ~ 1000 km (runaway growth) and then grow more slowly to ~ 3000- 
5000 km (oligarchic growth). Objects close to the central star grow faster. In roughly 1 Gyr, 
however, all objects reach a characteristic maximum radius of 3000-5000 km. 
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log Time (yr) 



Fig. 5. — As in Fig. H] for the slope of the cumulative size distribution. During runaway 
growth, the slope decreases rapidly from q > 10 to q ~ 4. As oligarchic growth proceeds, the 
slope asymptotically approaches q ~ 3. Noise in the plots results from counting statistics in 
the mass bins. 
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Fig. 6. — Variation of the median g 4 with r max for all annuli at all times in an ensemble of 
disks at 15-75 AU starting with most of the mass in the largest planetesimals and the f w 
fragmentation parameters. The legend indicates x m , the scaling factor for the initial surface 
density (eq. [4]). The typical inter-quartile range is roughly 5q iqr ~ 2 at log r max ~ 2.25-2.75 
and 5q iqr ~ 0.5-1.0 at log r max > 3. For all x m , the absolute lower bound in the slope is g 4 ~ 
2. The upper bound has g 4 roughly inversely proportional to r max and slowly decreasing 
with x m . In calculations with smaller x m , slower growth allows the slope to reach the lower 
bound when the largest objects are relatively small. 




3 3.5 
log Maximum Radius (km) 



Fig. 7. — As in Fig. [6] for calculations with larger planetesimals. The lower panel repeats 
Fig. |6j The legend indicates the initial radius of the largest planetesimals. Compared to 
calculations with 1 km planetesimals, calculations with 10-100 km planetesimals yield larger 
planets and larger g 4 at fixed log r max . 
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Fig. 8. — As in Fig. [7] for disks at 30-150 AU over a 10 Gyr evolution time. Results at 
30-150 AU yield the same relationship for g 4 as a function of r max and x m . Longer evolution 
times in calculations at 30-150 AU produce larger r max than the calculations at 15-75 AU. 
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Fig. 9. — As in Fig. [8] for disks at 30-150 AU starting with a mix of planetesimal sizes 
and the f w fragmentation parameters. Calculations with a broader initial size distribution 
yield less massive planets and a shallower relation between g 4 and r max . For a fixed r max , 
calculations with larger r have steeper size distributions. 
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Fig. 10. — As in Fig. [9] for calculations at 30-150 AU with the f w (lower panel), f s (mid- 
dle panel), and f vs (upper panel) fragmentation parameters. All calculations begin with a 
broad initial size distribution of planetesimals and r$ = 1 km. Calculations with weaker 
planetesimals yield less massive planets. 
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Fig. 11. — Variation of the median q\ with r max for calculations of disks at 30-150 AU, 
the f w fragmentation parameters, and r = 1 km (lower panel), 10 km (middle panel), and 
100 km (upper panel). Results for calculations with n c oc r -3 (upper branch of points in 
each panel) and n c oc r~ a5 (lower branch of points in each panel) are included in all panels. 
Colors indicate results for different values of Fig. [61 The evolution of the slope is 

remarkably independent of r . Until the onset of the collisional cascade at log r max w 2.5- 
3, calculations 'remember' the initial slope of the size distribution at 0.1-1 km. Once the 
collisional cascade begins, the slope rapidly evolves to q\ ~ 2. As r max reaches ~ 3000 km, 
the slope declines to q\ ~ 0.0-0.5. 
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Fig. 12. — As in Fig. [TT]for q 2 , the slope of the size distribution at 1-10 km. For calculations 
with small r , the evolution of qi is independent of the slope of the initial size distribution 
at smaller sizes; as r max grows from ~ 30 km to ~ 1000 km, q2 monotonically declines from 
q 2 5-6 to q2 ~ 4. When tq is larger, q 2 remains close to its initial value until the onset 
of the collisional cascade at log r max rs 3. When the cascade begins, q 2 evolves to ~ 2 and 
remains close to this value independent of r max . 



• •• - 

()l i I i I i L 




3 3.5 4 4.5 

log Maximum Radius (km) 



Fig. 13. — As in Fig. [IT] for g 3 , the slope of the size distribution at 10-100 km. For 
calculations with r < 10 km, the evolution of g 3 is independent of the slope of the initial size 
distribution at smaller sizes; as r max grows from ~ 300 km to ~ 1000 km, g 3 monotonically 
declines from g 3 « 5-6 to g 3 ~ 3-4. At any log r max < 3, the variation in g 3 with x m is 
larger for calculations with smaller tq. Once the collisional cascade begins, the slope slowly 
evolves to g 3 ~ 2-3, with g 3 w 2 favored at late times. 
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Fig. 14. — Fraction of calculations that achieve a target slope = 5 ± 0.25 and maximum 
size log r max = 2.80 ± 0.05 (left panels), log r max = 3.00 ± 0.05 (middle panels), or log r max 
= 3.15 ± 0.05 (right panels) as a function of the initial disk mass (x m ). In each panel, the 
legend indicates r . Colored lines show results for different starting conditions; violet: a = 
15-75 AU; n c oc r -0 - 5 , f w fragmentation parameters; blue: a = 30-150 AU; n c oc r~ - 5 , f w 
fragmentation parameters; green: a = 30-150 AU; n c oc r -3 , f w fragmentation parameters; 
orange: a = 30-150 AU; n c oc r~ 3 , f s fragmentation parameters; magenta: a = 30-150 AU; 
n c oc r -3 , f vs fragmentation parameters. For r = 1 km (lower panels), any range of a, most 
initial disk masses, and all fragmentation parameters, 10% to 30% of the calculations achieve 
the target slope and maximum radius; the success rate decreases slowly with r max . For r 
= 10 km (middle panels) and tq = 100 km (upper panels), 20% to 50% of calculations with 
larger disk masses (x m > 0.3) - but none of the calculations with smaller disk masses (x m < 
0.3) reach the target; the success rate increases slowly with r max . 
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Fig. 15. — As in Fig. [TUfor a target slope = 3 ± 0.25. For r = 1 km (lower panels), 
any range of a, most initial disk masses, and all fragmentation parameters, 20% to 50% of 
the calculations reach the target slope and maximum radius; the success rate increases with 
larger r max . For r = 10 km, calculations with large disk masses (x m > 0.3) reach the target; 
other calculations fail. For r = 100 km, none of the calculations match the target. In all 
calculations with r = 10-100 km, the success rate is independent of r max . 
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Fig. 16. — The median evolution time to reach r max = 1000 km as a function of semimajor 
axis for calculations at 15-150 AU with the f w fragmentation parameters. Left panels show 
results for calculations where all of the initial mass is mostly in large planetesimals (n c oc 
r -0,5 ). Right panels show results for calculations where the initial mass is equally distributed 
in log mass (n c oc r~ 3 ). In each panel, the legend indicates the initial radius of the largest 
planetesimal. Colored points indicate results for initial values of x m from Fig. |6j The 
black horizontal line illustrates the app r oximate time of the Late Heavy Bombardmen t (e.g., 



Tera et al.lll974j ; lHartmann et al.ll200d : IStoffler fc Rvderll200ll IChapman et all 120071 ). 



max 






X m> r 



any 





1000 km 









Fig. 17. — Schematic diagram for constraints on coagulation models. The left column 
lists target values for g 4 and r max . The right column lists the constraints on x m and r . 
Individually, measurements of g4 and r max for the cold (first row from top) or the hot 
(second row) population place few constraints on x m or r . Combined (third row), these 
data limit coagulation models to specific combinations of x m and r . When limits on the 
formation time are added (fourth row), successful models are limited to massive disks with 
r = 1 km. 



